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Abstract 

This  research  explores  an  innovative  sampling  method  nsed  to  condnct  nncer- 
tainty  analysis  on  a  system  with  one  random  inpnt.  Given  the  distribntion  of  the 
random  inpnt,  X,  we  seek  to  hnd  the  distribntion  of  the  ontpnt  random  variable  Y. 
When  the  fnnctional  form  of  the  transformation  Y  =  g{X)  is  not  explicitly  known, 
complicated  procednres,  snch  as  stochastic  projection  or  Monte  Garlo  simnlation, 
mnst  be  employed.  The  main  focns  of  this  research  is  determining  the  distribntion 
of  the  random  variable  Y  =  g{X)  where  g{X)  is  the  solntion  to  an  ordinary  dif¬ 
ferential  equation  and  X  is  a  random  parameter.  Here,  y  =  g{X)  is  approximated 
by  constrncting  a  sample  {Xj,  1/}  where  the  Xj  are  not  random,  bnt  chosen  to  be 
evenly  spaced  on  the  interval  [a,b]  and  Tj  =  g{Xi).  Using  this  data,  an  efficient 
approximation  ^(X)  5'(X)  is  constrncted.  Then  the  transformation  method,  in 

conjnnction  with  ^(X),  is  nsed  to  find  the  probability  density  fnnction  (pdf)  of  the 
random  variable  Y.  This  nniform  sampling  method  and  transformation  method  will 
be  compared  to  the  stochastic  projection  and  Monte  Garlo  methods  cnrrently  being 
nsed  in  nncertainty  analysis.  It  will  be  demonstrated,  throngh  several  examples,  that 
the  proposed  nniform  sampling  method  and  transformation  method  can  work  faster 
and  more  efficiently  than  the  methods  mentioned. 
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A  SAMPLING  AND  TRANSFORMATION  APPROACH  TO 


SOLVING  RANDOM  DIFFERENTIAL  EQUATIONS 

I.  Introduction 

1.1  Background 

Uncertainty  analysis  provides  us  with  information  regarding  the  uncertainty  of 
system  output  when  information,  about  inputs  into  the  system,  is  known.  If  we  know 
the  distribution  of  the  random  inputs,  X,  into  the  system,  then  we  seek  to  hud  the  dis¬ 
tribution  of  the  output  random  variable  Y.  The  complexity  of  this  procedure  depends 
on  whether  the  functional  form  of  the  transformation  Y  =  g{X)  is  known  explicitly 
or  not.  If  it  is  explicitly  known,  standard  statistical  techniques  can  be  applied.  If  it 
is  not  known,  then  other,  more  complicated  procedures  must  be  employed.  The  main 
focus  of  this  research  is  determining  the  distribution  of  the  random  variable  Y  =  g{X) 
where  g{.)  is  the  solution  to  an  ordinary  differential  equation  (ODE)  and  X  is  the  set 
of  random  parameters.  We  define  a  random  ODE  as  a  differential  equation  where  at 
least  one  parameter  or  initial  condition  is  a  random  variable.  Our  goal  is  to  introduce 
a  more  efficient  method  to  quantify  the  uncertainty  caused  by  random  inputs  into  a 
system.  This  method,  which  we  will  refer  to  as  the  uniform  sampling  method,  incor¬ 
porates  a  very  simple  approach  to  seemingly  complex  problems.  Here,  Y  =  g{X)  is 
approximated  by  constructing  a  sample  {Xj,  R}  where  the  Xi  are  not  random,  but 
chosen  to  be  evenly  spaced  on  a  set  determined  by  the  range  of  X  and  R  =  g{Xi). 
Using  this  data,  an  efficient  approximation  g{X)  k,  g{X)  is  constructed.  Then  an 
efficient  alternative  to  Monte  Carlo  simulation,  called  the  transformation  method,  is 
used  in  conjunction  with  g{X)  to  find  the  probability  density  function  (pdf)  of  the 
random  variable  Y . 

The  concept  of  random  input  into  a  system,  mainly  engineering  systems,  has 
been  studied  extensively  and  many  results  have  been  attained.  The  stochastic  projec- 
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tion  method,  described  later  in  this  text,  provides  a  way  to  determine  the  distribution 
of  the  output  random  variable  based  on  the  random  input.  To  build  an  expansion 
without  knowing  the  form  of  the  response,  a  random  basis  orthogonal  to  the  distribu¬ 
tion  of  the  input  uncertainty  is  typically  selected.  When  polynomials  are  selected  as 
the  basis,  the  method  is  referred  to  as  polynomial  chaos  [8].  The  concept  of  Homoge¬ 
neous  Chaos  was  pioneered  by  Norbert  Wiener,  in  1938,  involving  the  use  of  Hermite 
polynomials  and  the  input  random  variable  having  a  standard  normal  distribution  [4]. 
Throughout  the  text,  this  will  be  referred  to  as  Hermite-Chaos.  In  the  generalized 
polynomial  chaos,  called  the  Askey-Chaos,  the  polynomials  are  not  restricted  to  the 
Hermite  polynomials  and  the  random  variables  are  not  restricted  to  be  Gaussian  ran¬ 
dom  variables  [9].  When  using  the  polynomial  chaos  method,  it  is  standard  practice 
to  use  polynomials  that  are  orthogonal  with  respect  to  the  pdf  of  the  input  random 
variable.  In  Table  (1),  the  corresponding  type  of  polynomial  and  their  associated 
random  variables  are  listed  [9].  In  this  table,  the  support  set  is  the  region  where  the 
random  variable  is  dehned.  We  will  employ  Hermite-Chaos  in  this  research  since,  by 
invoking  the  central  limit  theorem,  the  Gaussian  distribution  appears  to  be  the  most 
likely  candidate  for  many  physical  applications  [8]. 


Random  Variables 

Orthogonal  Polynomials 

Support 

Continuous 

Gaussian 

Hermite 

(  — CXD,  OO) 

Gamma 

Laguerre 

[0,oo) 

Beta 

Jacobi 

[a,b] 

Uniform 

Legendre 

[a,b] 

Discrete 

Poisson 

Charlier 

{0,1,2,...} 

Binomial 

Krawtchouk 

(0,1,  2,. ..A) 

Negative  Binomial 

Meixner 

{0,1,2,...} 

Hypergeometric 

Hahn 

{0,1,  2,. ..A} 

Table  1:  Correspondence  of  Polynomials  and  Random  Variables  for  Dif¬ 
ferent  Askey-Chaos  (N  >  0  is  a  finite  integer) 
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While  Hermite-Chaos  is  useful  in  the  analysis  of  stochastic  processes,  efforts  have 
also  been  made  to  apply  it  to  model  uncertainty  in  physical  applications  [13].  These 
applications  include  structural  fatigue,  structural  reliability,  structural  mechanics, 
linear  structural  dynamics,  nonlinear  random  vibration,  soil  mechanics,  soil-structure 
interaction,  and  simulation  of  non-Gaussian  random  fields  [7].  For  example,  in  [9], 
Hermite-Chaos  is  used  to  determine  the  cross-flow  displacement  of  an  elastically- 
mounted  circular  cylinder  subject  to  stochastic  inputs.  The  cylinder  is  free  to  move 
in  the  y-direction  but  not  in  the  x-direction.  This  problem  has  two  random  inputs; 
damping  of  the  cylinder  and  the  stiffness.  In  this  research,  Hermite-Chaos  will  be  used 
on  linear  and  nonlinear  random  ordinary  differential  equations,  with  one  random 
input,  to  find  the  functional  form  of  the  transformation  Y  =  g{X)  when  X  is  a 
standard  normal  random  variable. 

Another  way  to  approximate  the  distribution  oiY  =  g{X)  is  through  the  use 
of  Monte  Carlo  simulation.  Here,  a  random  sample  from  the  input  random  variable 
is  used  to  approximate  the  distribution  of  the  output  random  variable.  This  can  be 
slow  with  respect  to  computational  time,  but  produces  accurate  results.  In  Xiu  and 
Karniadakis  [4],  it  was  shown  that  Hermite-Chaos  was  substantially  faster  than  Monte 
Carlo  simulations  for  low  dimensional  stochastic  inputs. 

Modelling  and  simulation  complexity  challenges  the  state  of  computing  and  will 
continue  to  demand  improvements  in  accuracy  and  efficiency  [12].  Many  of  the  ma¬ 
jor  modelling  and  simulations  efforts,  supported  by  the  U.S.  Department  of  Energy 
(DOE),  must  address  uncertainty  and  reliability  in  order  to  provide  confidence  levels 
for  the  results  [12].  For  example,  reliability  in  predictions  is  crucial  because  ulti¬ 
mately  these  predictions  can  impact  risk  assessment  decisions  [12].  Predictive  models 
of  microbial  cell  and  community  functions  are  hindered  by  the  lack  of  reliable  and 
complete  data  for  estimating  kinetic  parameters  and  identifying  cell  signaling  path¬ 
way  structure  [12].  Answering  this  question  is  a  key  to  the  success  of  Genomes  to  Life 
as  the  results  of  combined  experimental  and  computational  studies  on  microbial  sys- 
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terns  are  used  to  design  strategies  for  bioremediation  of  the  nuclear  weapons  complex, 
converting  biomass  to  fuels,  and  carbon  sequestration  [12], 

In  a  proposal  submitted  to  the  DOE  from  Pacihc  Northwest  National  Labo¬ 
ratory,  a  method  referred  to  as  Smart  Sampling  was  introduced.  Smart  sampling 
implies  the  choice  of  an  experimental  design  to  specify  a  selection  of  input  parame¬ 
ters  that  will  best  cover  the  possible  outputs  of  the  simulation  [12].  A  number  of  these 
sampling  methodologies  have  been  drawn  from  the  design  of  physical  experiments  and 
modihed  for  the  design  of  computer  experiments  (This  modihcation  is  needed  because 
the  notions  of  blocking,  replication,  and  randomization  do  not  apply  to  deterministic 
computer  experiments)  [12].  This  smart  sampling  approach  was  used  to  determine 
efficient  input  into  a  computationally  intensive  computer  model.  The  computer  model 
can  be  thought  of  as  a  ’’black  box”  which  can  be  evaluated  but  nothing  else  is  known 
about  the  system.  Then,  by  the  sampling  method,  a  functional  approximation  to  the 
’’black  box”  can  be  constructed  and  used  to  determine  the  distribution  of  the  out¬ 
put.  The  overall  impact  of  the  project  is  the  enabling  of  accurate  approximation  of 
mean  behavior  and  uncertainty  in  complex  computational  models  at  a  greatly  reduced 
computational  cost  relative  to  Monte  Carlo  simulation  [12]. 

In  a  recent  research  project  involving  a  pitch  and  plunge  airfoil  with  cubic 
and  pentic  structural  parameters  and  uncertainties  in  the  spring  constant  and  initial 
pitch,  a  very  efficient  algorithm  for  determining  the  propagation  of  uncertainties  in 
a  highly  nonlinear  system  was  presented  [6].  These  uncertainties  were  allowed  to 
propagate  in  a  time  dependent  manner  until  a  limit  cycle  oscillation  or  a  dynamically 
stable  solution  was  achieved  [6].  The  stochastic  algorithm  consisted  of  building  an 
interpolating  function  in  the  stochastic  domain  through  sampling  and  the  use  of  a 
multivariate  B-spline  [6].  Advantages  of  this  method  were  that  expected  values  were 
never  computed  or  stored  and  two  orders  of  magnitude  efficiency  in  computational 
time  over  a  traditional  Monte  Carlo  simulation  were  obtained  [6].  Based  on  the 
PDFs  predicted  by  the  stochastic  algorithm,  a  very  rapid  and  accurate  estimate  of 
the  probability  of  failure  was  obtained  [6]. 
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1.2  Problem  Definition 

While  Hermite-Chaos  and  Monte  Carlo  simnlation  are  effective  methods  of  de¬ 
termining  the  distribntion  of  the  ontpnt  random  variable  of  some  systems,  there  are 
problems  when  neither  one  prodnces  satisfactory  resnlts.  For  complex  models  that 
consnme  considerable  computational  resources,  the  Monte  Carlo  method  may  be  very 
costly  and  inefficient.  Also,  this  approach  does  not  readily  provide  information  about 
the  sensitivity  of  model  outputs  to  specific  parametric  uncertainties  [10].  In  addition, 
polynomial  chaos  can  fail  when  the  transformation  g  :  X  ^  Y  is  not  smooth.  This 
occurs  because  polynomial  chaos  approximates  the  transformation  using  polynomial 
basis  functions  with  global  support.  Classic  results  from  polynomial  approximation 
theory  indicate  that  in  order  to  approximate  discontinuities  well,  piecewise  polyno¬ 
mials  must  be  used  [11]. 

In  the  following  discussion,  the  three  methods,  Hermite-Chaos,  Monte  Carlo 
simulation,  and  the  uniform  sampling  method  will  be  demonstrated  and  compared 
through  the  solution  of  both  linear  and  nonlinear  random  ODEs.  Computation  time 
and  comparison  to  the  actual  transformation  will  be  analyzed  using  these  three  meth¬ 
ods.  Therefore,  the  validity  of  the  uniform  sampling  method  will  be  determined  based 
on  the  results  of  these  tests. 
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II.  Methodology 

There  are  several  possible  approaches  for  finding  the  distribution  of  a  random  vari¬ 
able  y,  which  depends  on  a  random  variable  X  when  the  distribution  of  X  is  known. 
These  approaches  can  be  divided  into  two  cases.  The  first  case  is  when  the  functional 
form  of  the  transformation  g  :  X  ^  Y  is  known  explicitly.  The  second  case  is  when 
the  functional  form  of  g{.)  is  not  known  explicitly,  but  g{.)  can  be  evaluated  or  at 
least  approximated.  This  includes  the  situation  when  g  is  the  solution  to  a  differen¬ 
tial  equation  but  an  explicit  expression  for  g  is  unknown.  In  this  case,  a  numerical 
method  could  be  used  to  approximate  g.  These  ideas  will  be  explained  in  the  next 
several  sections  and  are  summarized  in  Table  (2). 


Approximate  g(X) 

Find  Distribution  of  Y 

Functional  Form  of  g(X) 
is  Known  Explicitly 

N/A 

Transformation  Method  or 
Monte  Carlo  Simulation 

Functional  Form  of  g(X) 
is  Not  Known  Explicitly 

Stochastic  Projection  or 
Uniform  Sampling  and  Interpolation 

Transformation  Method  or 
Monte  Carlo  Simulation 

Table  2:  Summary  of  Methods  Used  to  Find  Transformations  and  Distribn- 
tions 


2.1  The  Functional  Form  of  g{X)  is  Known  Explicitly 

In  this  section,  we  will  use  the  transformation  method  and  Monte  Carlo  Simula¬ 
tion  to  determine  the  pdf  of  the  system  output  random  variable  Y.  For  the  following 
examples  Y  =  sin^{X)  will  be  used  as  the  known  transformation  g{X). 

2.1.1  Finding  the  PDF  With  the  Transformation  Method.  Let 

X  be  a  random  variable  with  known  probability  density  function  (pdf)  fx{x)  with 
support  X  and  cumulative  distribution  function  (cdf)  Fx{x)  =  P{X  <  x).  Also,  let 
Y  =  g{X)  be  a  random  variable  that  depends  on  X.  If  g  and  the  pdf  of  X  are  known. 
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then  the  cdf  of  Y,  Fyiu),  can  be  determined  as  follows. 


Fy{y)  =  P{Y  <y)  =  P{g{X)  <  y) 

=  P{{x  e  df  :  g{x)  <  y]) 

=  /  fx{x)  dx 

J  {x&X:g{x)<y} 

The  pdf  of  Y  is  fy{y)  =  ^Fy{y).  This  approach  for  determining  the  distribution  of 
Y  will  be  referred  to  as  the  transformation  method.  As  discussed  in  [3]  this  method 
requires  g  to  be  monotone  over  a  partition  of  A. 

The  following  example  from  Casella  and  Berger  [3]  illustrates  this  method.  Suppose 
X  has  a  uniform  distribution  on  (0,27r), 


fx{x) 


A  0  <  X  <  27r 

Ztt 

0  otherwise 


and  let  g{x)  =  sin^(x)  or 


Y 


sin^(X)  0  <  X  <  27r 
0  otherwise 


Y=SinA 


Figure  1:  Transformation  Y  =  sin^(X)  Where  X  is  Uniform(0,27r) 
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Then,  referring  to  Figure  (1) 


Fy^y)  =  P(Y  <  y)  =P{X  <  Xi)  +  P(X2  <  X  <  X3)  +  P(X4  <  X) 
=  ■■■  =  2  [P{X  <  xi)  +  P(X2  <  X  <  vr)] 


where  xi  <  X2  and  xi,X2  are  the  solutions  to  sin^(x)  =  y,  0  <  x  <  ti.  Then  xi  is  the 
smallest  x  >  0  such  that  x  =  arcsin(y/^)  and  X2  =  vr  —  Xi.  So 


My) 


Xi  TT  —  X2 

^  ^  271 


2xi 

71 


2 

71 


arcsin(-y/y) 


and 


fviy) 


d  [2 
dy  [tt 


arcsin(y^) 


1 

T^\/yO--y) 


In  Figure  (2)  the  pdf  and  cdf  of  X  are  shown.  In  Figure  (3)  the  pdf  and  cdf  of 
Y  =  sin^(X)  are  shown. 

fx  Fx 


Figure  2:  PDF  and  CDF  of  X  uniform(0,27r) 


This  method  relies  on  being  able  to  partition  X  into  regions  on  which  g{X)  is 
monotone.  It  also  requires  that  for  a  given  Y,  X  =  g~^{Y)  can  be  computed.  If 
either  of  these  requirements  are  not  met  or  is  especially  expensive  to  compute,  this 
method  will  fail. 


Figure  3:  PDF  and  CDF  of  F  =  sin^(X)  with  X  uniform(0,27r) 


y 

Figure  4:  Histogram  of  F  =  sin^(X)  with  X  uniform(0,27r) 

2.1.2  Finding  the  PDF  With  the  Monte  Carlo  Method.  The  trans¬ 
formation  y  =  g{x)  can  be  evaluated  for  any  x.  Also  assume  the  distribution  of  X  is 
known  or  at  least  a  method  of  generating  a  random  sample  from  X  is  available.  One 
way  to  approximate  the  pdf  of  F  =  g{X)  is  to  generate  a  random  sample  {xi, . . . , 
from  X  and  then  substitute  the  sampled  values  into  y  =  g{x),  producing  a  random 
sample  {yi, . . .  ,?/„,}  from  F.  A  histogram  of  these  sample  values  approximates  the 
pdf  of  F.  The  result  of  doing  this  with  g[x)  =  sin^  x  and  X  uniform(0,  27r)  is  shown 
in  Figure  (4),  which  matches  the  pdf  in  Figure  (3)  quite  well.  There  are  several  other 
ways  to  approximate  the  pdf  of  F  given  a  random  sample  {yi, . . . ,  yn}  which  are  out- 
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lined  in  [2],  Of  course,  these  approaches  can  be  applied  to  any  sample  {yi, . . . ,  yn} 
whether  it  was  generated  as  y  =  g{x)  or  not.  These  methods  will  be  referred  to  as 
Monte  Carlo  methods  because  they  are  based  on  a  random  sample.  While  Monte 
Carlo  simulation  is  the  standard  method  in  practice  today,  it  has  some  limitations.  It 
may  be  costly  and  inefficient  especially  when  the  system,  under  analysis,  is  compu¬ 
tationally  expensive.  To  obtain  a  useful  pdf,  the  histogram  must  be  smoothed.  The 
transformation  method,  however,  produces  a  more  accurate  and  usable  pdf  without 
any  smoothing  or  extra  steps. 

2.2  The  Functional  Form  of  g{X)  is  Not  Known  Explicitly 

In  this  section,  we  will  use  the  transformation  method  and  Monte  Carlo  Simula¬ 
tion  to  determine  the  pdf  of  the  system  output  random  variable  Y.  For  the  following 
examples  Y  =  sin^{X)  will  be  used  as  the  unknown  transformation  g{X).  That  is, 
the  uniform  sampling  method  and  Hermite-chaos  will  be  used  to  approximate  the 
transformation  Y  =  sin^{X). 

In  Table  (2),  there  are  two  methods  to  obtain  the  approximation  to  the  trans¬ 
formation  g{X),  stochastic  projection  and  uniform  sampling.  To  demonstrate  the 
uniform  sampling  method,  we  suppose  g{X)  =  sin^{X)  where  X  has  a  uniform  dis¬ 
tribution  on  (0,  27r).  This  transformation  can  be  approximated  by  uniformly  sampling 
X  over  the  interval  (0,  27r)  producing  a  sample  {xi, . . .  Then  let  j/j  =  g{xi),  pro¬ 

ducing  a  collection  of  points  {a:*,  ?/*}  where  i  =  1, . . . ,  n.  Then  we  can  interpolate  g  at 
{xi^yi}  producing  the  approximation  g{X)  ^  g{X).  Figure  (5)  shows  the  piecewise 
linear  interpolation  of  sin^{X)  using  100  evenly  spaced  samples. 

Next,  we  demonstrate  the  stochastic  projection  method  for  the  transformation 
g{X)  =  sin^{X)  where  X  has  a  Uniform  distribution  on  (0,27r).  From  Table  (1), 
we  see  that  the  Legendre  polynomial  is  the  correct  orthogonal  polynomial  to  use 
with  random  variables  from  the  Uniform  distribution.  The  result  in  Figure  (6)  shows 
that  the  stochastic  projection  method  works  quite  well.  The  actual  transformation, 
g{X)  =  sin^(X),  is  shown  as  the  solid  plot  while  the  approximated  transformation. 
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ghat 


Figure  5:  g{X)  Using  Uniform  Sampling  and  Interpolation  With  100  Samples 

g{X),  is  the  dashed  plot.  The  stochastic  projection  method  will  be  discussed  in  the 
next  chapter. 

ghat 


Figure  6:  Approximate  and  Actual  Transformation  of  U  =  g{X)  Using  Polynomial 
Chaos  of  Degree  8  With  X  Uniform(0,27r) 


2.2.1  Finding  the  PDF  With  The  Transformation  Method.  Once 
g{X)  is  obtained,  the  transformation  method  can  be  applied  to  get  the  cdf  and  pdf  of 
the  output  random  variable.  The  results  for  g{X)  produced  by  sampling  are  found  in 
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Figure  (7).  Notice  that  the  approximation  to  g{X)  using  sampling  produces  almost 
the  exact  pdf  and  cdf  of  the  output  random  variable.  The  error  plot  shown  in  Fig¬ 
ure  (8),  calculated  as  \g{X)  —  g{X)\,  verifies  that  the  g{X)  produced  via  sampling  is 
a  good  approximation  to  the  actual  transformation.  Similar  results  can  be  achieved 
using  g{X)  produced  by  stochastic  projection. 

Fy  fy 


Figure  7:  PDF  and  CDF  oi  Y  =  g{X)  Using  Transformation  Method  with  X 
uniform(0,27r) 


Error 


Figure  8:  Error  Plot  Y  =  \pdf  —  pdf\  Using  Transformation  Method  with  X 
uniform(0,27r) 
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2.2.2  Finding  the  PDF  With  the  Monte  Carlo  Method.  Using  the 
Monte  Carlo  method  described  in  section  2.1.2,  we  can  approximate  the  pdf  of  the 
output  random  variable  using  g{X),  found  in  section  2.2,  with  X  having  a  uniform 
distribution  on  (0,27r).  Once  the  sample  is  generated,  a  histogram  is  constructed 
to  form  the  pdf.  The  pdf  is  shown  in  Figure  (9)  where  g{X)  was  generated  by 
the  uniform  sampling  method  with  100  samples.  As  expected,  this  function  mirrors 
the  exact  pdf  in  section  (2.1).  Also,  if  evaluation  of  the  function  g{x)  is  expensive, 
generating  a  sample  {yi, . . . ,  ?/„}  from  {xi, . . . ,  using  Monte  Carlo  will  also  be 
expensive.  Since  g  is  typically  a  piecewise  polynomial  and  can  be  evaluated  efficiently, 
one  approach  would  be  to  generate  a  random  sample  {xi, . . .  ,x„}  using  Monte  Carlo 
and  set  iji  =  g{xi)  to  efficiently  create  the  sample  {^j, . . .  ,^n}-  Then  a  histogram 
of  this  sample  approximates  the  pdf  of  Y .  Also,  since  the  pdf  of  X  is  known,  the 
transformation  method  using  g  could  be  used  to  better  approximate  the  cdf  and  pdf 
of  U. 

fy 

8 

6 


Figure  9:  PDF  of  Y  =  g{X)  Using  Monte  Carlo  With  X  uniform(0,27r) 
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III.  Random  Ordinary  Differential  Equations 

The  polynomial  chaos  and  nniform  sampling  methods  can  be  applied  directly  to  ran¬ 
dom  ordinary  differential  eqnations  (ODEs).  In  this  case,  the  explicit  fnnctional  form 
of  the  transformation  is  not  known,  bnt  certain  characteristics  of  the  system,  snch 
as  a  rate  of  change,  are  known.  Generally,  systems  nnder  analysis  exhibit  random 
inpnts.  A  random  ODE  is  one  where  the  differential  operator  L[u{t)]a\  depends  on 
a  parameter  set  a  and  one  or  more  elements  of  a  are  random  variables.  In  addition, 
the  initial  or  bonndary  conditions  conld  be  random  variables.  The  methods  nsed  to 
approximate  the  transformation  will  be  illnstrated  with  two  examples,  a  linear  ODE 
and  a  nonlinear  ODE. 

3.1  Linear  Example 

The  hrst  example  is  the  linear  ODE 

mu  +  cit  +  ku  =  0 

n(0)  =  0  (1) 

m'(0)  =  V 

which  models  the  spring  mass  damper  system  fonnd  in  Fignre  (10).  Here  k  is  the 
spring  constant,  m  is  the  mass,  and  c  is  the  damping  coefficient.  The  solntion  to  (1) 
is  emphasizing  that  u  depends  on  the  independent  variable  and  five 

parameters.  This  can  be  viewed  as  a  transformation  Y{t)  =  g{X).  If  any  com¬ 
ponent  of  X  is  random,  then,  for  hxed  t,  Y  (t)  is  a  random  variable  and  the  dis- 
tribntion  of  Y{t)  can  be  determined  nsing  the  methods  of  Chapter  2  provided  the 
distribntion  of  the  random  components  of  X  are  known.  Here,  for  u{t-,X),  we  have 
X  =  {m,c,k,u{0),u\0)). 

Of  the  methods  presented  in  Chapter  2  for  determining  the  distribntion  of  Y [t), 
we  will  focns  on  the  case  where  the  fnnctional  form  of  g{.)  is  not  known  explicitly. 
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Figure  10:  Spring  Mass  Damper  System 


Even  though  the  model  problem  (1)  can  be  solved  explicitly,  explicit  solutions  are  not 
available  for  many  problems  of  interest  in  applications. 

3.1.1  Approximating  the  Transformation  g-.X^Y  Using  the  Stochas¬ 
tic  Projection  Method.  We  begin  with  stochastic  projection  which  is  described 
extensively  in  Xiu  and  Karniadakis  [4]  and  in  Ghanem  and  Spanos  [8].  The  approxi¬ 
mate  transformation  is  represented  by 


p 


(2) 


Given  0*,  we  will  show  how  to  find  Xi.  How  the  fnnctions  (ji  are  chosen  will  be 
discussed  later  in  this  chapter.  The  first  step  is  to  find  xft)  for  i  =  0, . . .  ,P  where 
P  is  the  order  of  the  polynomial  expansion.  Given  the  model  problem  in  Eqnation 
(1),  we  consider  the  case  where  there  is  one  random  parameter,  which  is  the  spring 
constant  k.  We  let  k  be  defined  by 


p 


(3) 
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Next,  we  substitute  Equation  (3)  and  (2)  into  (1).  This  results  in 


p  p  ^  p  p 

—  =  0.  (4) 

2  =  0  2=0  2=0  j  =  0 

Then  we  multiply  equation  (4)  by  (pi  where  I  =  0, . . . ,  P  and  a  weight  function  w{x) 
and  integrate  on  [a,b].  Letting  <  f,g  >  =  la  faw  dx,  we  get 

<(ph(pi>  +  ^  EJIo  <  (pi,  (l^i>  +  ^  E!1o  Ef=o  kiXj(t)  <  (pi(pj,  (pi>  =  0 

for  /  =  0, . . . ,  P. 

(5) 

If  we  choose  (p  so  they  are  orthogonal  with  respect  to  <  .  >,  then  the  result  is 


[p(0  +  <  00  E^lo  Ef=0  kiXj{t)  <  (piCpj,  ^i>  =  0 

for  /  =  0, . . . ,  P. 


Dividing  through  hy  <  (pi,  (pi  >,  we  get 


P(0  +  ^Xl(t)  +  EL  Ef=0  kiXjit)  <  (Pi(Pj,  (Pi>  =  0 

for  Z  =  0, . . . ,  P. 

To  convert  this  system  of  second  order  ODEs  to  a  system  of  first  order  ODEs,  we  let 


Ui  =  Xi(t) 

ill  =  xi(t)  =  vi 

vi  =  xiit)  =  -^Vi-  Ef=o  Ef=o  <  0*0T  0;  > 

for  /  =  0, . . . ,  P. 
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This  gives  us 


Mo 

^^0 

Ml 

Vi 

M  p 

,  b  = 

Vp 

Vo 

m'^0  m<</.0,</>0>  ^*=0  ^i=0  <  4>i4’jy  00  > 

Vi 

vp 

m<4>p,4,p>  Ylii=oYlij=0^iVj  <  4>i4>j',4>P  > 

The  resulting  system  that  needs  to  be  solved  is 


z  =  Sz 


(10) 


where 


S  = 


1  .00 

0  .00 


0 


0 


m<(f>o,(t>o>  Ei=0  ^  4’0  >)  ■••)  X]i=0  ^  4>i4>Pi  4>0  >] 


m<ipp,<l>p>  Ei=0  ^  4*i4'Oi  bp  >I  ■■■^  Si=0  ^  4>i4‘P^  4>P  >] 


_  ^ 
m 

0 

0 

0 

0 


1  0 
0  1 
0  0 
0 
0 
0 
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To  simplify  the  system,  we  choose  the  (f)j  such  that  the  inner  product,  <  4>j  >  =  1 

for  all  j  =  0, . . . ,  P.  We  also  set  each  ki  =  0  for  i  >2.  Then  the  matrix  S  will  be 


(12) 


Thus,  the  following  description  can  be  used  to  solve  the  system  of  ODEs  to  obtain 
the  approximate  transformation  using  stochastic  projection.  We  define  the  following: 

•  a,b  -  the  interval  on  which  the  random  variable  ^  is  defined 

•  (j)i{^)  -  stochastic  basis  function 

•  m,c  -  parameters  in  ODE 

•  k  -  spring  constant 

•  ko  -  mean  value  of  k  (will  be  denoted  k) 

•  ki  -  standard  deviation  of  k  (will  be  denoted  k) 

•  V  -  initial  velocity 

•  tf  -  hnal  time 

•  w{^)  -  the  pdf  of  the  distribution  of  the  input  random  variable 

We  must  construct  the  matrices  A  and  B  in  the  ODE  system 


X  -\r  cx  +  ( — /  H - A  ^B)x  =  0 


(13) 


m  m 


where  the  matrix  A  =  [<  >]•  These  inner  products  are  calculated  by 


(14) 


18 


The  matrix  B  is  given  by  B  =<  >  which  is  similar  to  the  A  matrix. 

These  inner  prodncts  are  calculated  by 


Once  we  have  A  and  B,  We  use  them  to  solve  the  system 


(15) 


z  =  Sz 


(16) 


where  S  is  the  matrix  defined  by 


S  = 


(17) 


The  result  of  solving  Equation  (16)  gives  us  Xk{t).  The  transformation  can  now  be 
approximated  using  the  equation  =  J2k=o^kit)(f>k{0- 


We  now  look  at  the  following  example  of  the  linear  ODE  described  in  (10).  For 
this  example,  the  polynomial  chaos  method  described  uses  the  one-dimensional  Her- 
mite  polynomials  in  terms  of  standard,  normally  distributed  random  variables  (.^). 
We  begin  by  defining  as  the  Hermite  polynomial  of  degree  i  and  w{^)  to  be  the 
pdf  of  the  standard  normal  distribution.  We  also  use  the  following  parameters  for  the 
model  problem  (1) 


a 

b 

m 

k 

k 

c 

V 

t 

-6 

6 

10 

10 

2 

4 

25 

10 

Table  3:  Parameters  for  Linear  Model  ODE  (1) 
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For  the  function  we  used  the  scaled  Hermite  polynomial 


MO  = 


(18) 


The  scaling  is  needed  so  that  <  (t)i,c()i  >=  1.  For  the  function  wM,  we  used 


w{0  = 


1 

e  2 


When  P  =  8,  the  matrices  A  and  B  are  as  follows 


(19) 


A  = 


ioooooooo\ 
010000000 
001000000 
000100000 
000010000 
000001000 
000000100 
000000010 
000000001 


(20) 


B  = 


0  1 

1  0 

0  1.41421 

0  0 


0 

1.41421 

0 

1.73205 

0 

0 

0 

0 

0 


0 

0 

1.73205 

0 

2 

0 

0 

0 

0 


0 

0 

0 

2 

0 

2.23607 

0 

0 

0 


0 

0 

0 

0 

2.23607 

0 

2.44949 

0 

0 


0 

0 

0 

0 

0 

2.44949 

0 

2.64575 

0 


0 

0 

0 

0 

0 

0 

2.64575 

0 

2.82843 


0 

0 

0 

0 

0 

0 

0 

2.82843 

0 


(21) 
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The  resultant  approximate  and  actual  transformations  are  given  in  Figure  (11). 
The  actual  transformation  was  found  using  Laplace  Transforms  on  the  differential 
equation 


X  + 


c  k 

— X  H - X  =  0. 

m  m 


The  exact  transformation  is 


(22) 


2ve  2m 

— .  sink 

\le-Hk  +  ik) 


—  {k  +  ^k) 
m 


(23) 


The  error  plot  shown  in  Figure  (12),  calculated  as  \g{0  ~  .^(01;  verifies  that 
the  g{^)  produced  with  polynomial  chaos  of  degree  8  is  a  good  approximation  to  the 
actual  transformation.  The  error  decreases  as  the  degree,  of  the  Hermite  polynomial 
used  in  the  expansion,  increases.  The  Mathematica  code  used  for  the  Hermite-chaos 
method  is  found  in  Appendix  A.l. 


g 


Figure  11;  Approximate  (dashed)  and  Actual  (solid)  Transformation  oiY  =  g{^) 
Using  Hermite-chaos  of  Degree  8. 


While  Hermite-Chaos  was  shown  to  be  an  effective  method  of  approximating 
the  transformation  for  the  linear  and  nonlinear  ODE  examples,  polynomial  chaos, 
in  general,  can  fail  when  the  transformation  g  :  X  —y  Y  is  not  smooth.  For  ex- 
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Error 


0.4 


0.3 


6 


-4  -2 


2 


4 


Figure  12:  Error  Plot  Y  =  —  ^(01  Using  Hermite-chaos  of  Degree  8. 

ample,  the  transformation  g{X)  =  tan{X),  where  X  is  a  random  variable  from  the 
Uniform{0,  n)  distribution,  is  discontinuous  at  a;  =  |.  When  using  polynomial  chaos, 
with  Legendre  polynomials  as  the  basis  functions  determined  from  Table  (1),  the  re¬ 
sults,  in  Figure  (13),  are  what  we  expect,  but  do  not  approximate  the  transformation 
well.  By  using  the  uniform  sampling  method,  described  in  the  next  section,  for 
g{X)  =  tan{X),  Figure  (13)  shows  that  the  uniform  sampling  method  approximates 
the  true  transformation  quite  well. 

ghat  ghat 


-1 


1 


Figure  13:  Approximate  (dashed)  and  Actual  (solid)  Transformation  oiY  =  tan{X) 
Using  Polynomial  Chaos  of  Degree  8  and  Uniform  Sampling  Method  With  80  Samples 
Where  X  Uniform(0,7r) 
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3.1.2  Approximating  the  Transformation  g-.X^Y  Using  the  Uni¬ 
form  Sampling  Method.  The  second  method  to  approximate  the  transforma¬ 
tion  y  =  g{X)  is  by  uniformly  sampling  the  transformation  and  using  interpolation. 
Approximate  the  transformation  g{x)  using  a  set  {{xi,  g{xi)}.  In  this  case,  the  set 
{xi, . . .  ,Xn}  is  not  random  but  is  chosen  to  be  evenly  spaced  on  the  interval  [a,b]. 
Then,  {g{xi), . . .  ,g{xn)}  is  calculated  by  solving  the  ode  system  for  each  Xi.  Using 
this  data,  an  approximation  g{x)  ~  g{x)  is  constructed.  The  approximated  and  ac¬ 
tual  transformation  are  shown  in  Figure  (14)  where  g{X)  is  the  solid  plot  and  g{X)  is 
the  dashed  plot.  The  Mathematica  code  used  for  the  uniform  sampling  method  can 
be  found  in  Appendix  A.2. 


g 


Error 


Figure  14:  Actual  and  approximate  Transformation  Y  =  g{X)  Using  Uniform  Sam¬ 
pling  Method  With  500  Samples  and  Error  Plot  Y  =  \g{X)  —  g{X)\ 


3.1.3  Finding  the  Distribution  With  Approximate  Transformation. 

Once  the  approximation  of  the  transformation  y  =  g{X)  is  found,  the  distribution 
of  the  output  variable  can  be  obtained.  The  methods  employed  to  accomplish  this  are 
Monte  Carlo  and  the  transformation  method  both  described  previously  in  Table  (2). 

3. 1.3.1  Monte  Carlo  Method.  The  distribution  of  the  output  vari¬ 
able  can  be  determined  using  the  ’’input”  Monte  Carlo  method.  This  is  a  technique 


23 


for  approximating  the  pdf  using  a  histogram  generated  by  producing  random  inputs 
and  solving  the  system  for  each  random  input.  The  result  is  obtained  in  Figure  (15). 


Figure  15:  PDF  of  Y  Using  ’’Input”  Monte  Carlo  With  20,000  Runs 


The  distribution  of  the  output  variable  can  also  be  determined  using  the  ’’out¬ 
put”  Monte  Carlo  method.  This  is  a  technique  for  approximating  the  pdf  by  approx¬ 
imating  the  transformation  and  use  Monte  Carlo  with  g{X)  evaluating  the  random 
inputs.  The  result  is  shown  in  Figure  (16).  The  Mathematica  code  for  Monte  Carlo 
simulation  can  be  found  in  Appendix  D. 

Fy 


Figure  16:  PDF  of  Y  Using  ’’Output”  Monte  Carlo  With  20,000  Runs  Where  g{x) 
Obtained  by  Hermite-Chaos  of  Degree  8. 
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3. 1.3. 2  Transformation  Method.  Using  g{x),  the  distribution  of 
the  output  variable  can  also  be  determined  using  the  transformation  method.  First, 
a  sample  X  =  {xi, . . .  ,x„}  is  generated.  This  sample  is  random  for  the  stochastic 
projection  method  and  equally  spaced  over  the  interval  [a,  b]  for  the  uniform  sampling 
method.  Then,  we  evaluate  Y  =  g{x)  =  {yi, . . .  ,yn}  producing  the  sample  Xi,yi. 
Next,  we  go  through  the  list  of  increasing  y  values,  hnding  the  limits  of  integration 
and  integrating  to  generate  the  data  to  build  the  cumulative  distribution  function 
(cdf).  Then  the  cdf  is  differentiated  to  get  the  pdf.  Both  the  cdf  and  pdf  are  shown 
in  Figure  (17).  The  Mathematica  code  for  the  transformation  method  can  be  found 
in  Appendix  C. 

Fy  fy 


Figure  17:  CDF  and  PDF  of  Y  Using  Transformation  Method  With  g{X)  Obtained 
by  Uniform  Sampling  With  512  Samples 

The  methods  used  to  obtain  the  distribution  of  the  output  random  variable  of  a 
system  with  one  random  input  are  outlined  in  Table  (2).  The  approximated  transfor¬ 
mation  g{X)  was  found  for  the  linear  ODE  with  the  uniform  sampling  method  and 
polynomial  chaos.  In  chapter  4,  the  accuracy  and  efficiency  of  the  uniform  sampling 
method  and  the  Hermite-chaos  method  will  be  analyzed  by  comparing  the  resulting 
approximate  transformations  to  the  true  transformation.  In  addition,  we  will  use 
Monte  Carlo  simulation  and  the  transformation  method  to  obtain  the  pdfs  using  the 
approximated  transformations  determined  from  Hermite-chaos  and  uniform  sampling. 
These  pdfs  will  also  be  compared  to  the  actual  pdfs  to  determine  the  accuracy  and 
efficiency  of  these  methods. 
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3. 2  Nonlinear  Example 

In  this  section,  we  consider  a  nonlinear  problem.  The  equation  of  the  spring 
mass  system  is 


mil  +  kill  +  k2U  +  (fcs  +  ^k^)u^  =  0 

m(0)  =  w  (24) 

m'(0)  =  V 

which  is  similar  to  the  linear  case  except  now  there  is  a  nonlinear  cubic  spring  stiffness. 
Since  g{X)  is  not  known  in  this  case,  it  must  be  approximated.  The  approximation  to 
the  transformation  will  be  obtained  by  the  stochastic  projection  method  and  the  uni¬ 
form  sampling  method.  Then,  the  transformation  method  or  Monte  Carlo  Simulation 
will  be  used  to  obtain  the  distribution  of  the  output  random  variable.  The  following 
computations  were  carried  out  with  ^  being  a  standard  normal  random  variable  and 
the  following  parameters  for  the  model  problem  (24).  The  actual  transformation, 
that  will  be  used  for  comparison,  was  obtained  by  sampling  using  8192  samples.  The 
result  is  shown  in  Figure  (18). 


t 

m 

ki 

k2 

h 

h 

w 

V 

10 

10 

0 

10 

1 

0.02 

10 

0 

Table  4:  Parameters  for  Nonlinear  Model  ODE  (24) 
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g 


Figure  18:  Actual  Transformation  by  Uniform  Sampling  Method  Using  8192  Sam¬ 
ples 


3.2.1  Approximating  the  Transformation  g-.X^Y  Using  the  Stochas¬ 
tic  Projection  Method.  This  method  follows  the  same  idea  given  for  the  linear 
ODE  described  previously.  The  method  is  shown  by  the  following  formulations. 

Consider  the  model  problem  from  (24).  The  approximate  transformation  is 
given  by 


p 

u{t,0  = 

i=0 


(25) 


and  we  dehne  k  to  be 


k  =  ks  +  where  ^  ~  A^(0, 1). 


(26) 


Next,  we  substitute  Equations  (26)  and  (25)  into  (24).  This  results  in 

P  P  _  ~  p 

'^Xi{t)(j)i  +  —  '^Xi{t)(j)i  +  (^3  +  ^h)  =  Q_  (^27) 

m  m 

Z=0  2=0  2=0 

Then  we  multiply  by  (f)j,  where  j  =  0, . . . ,  P,  and  a  weight  function  w{^)  and 
integrate  on  [a,  6].  Letting  <  f,g>  =  fiOgiO^iO  we  get 
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(28) 


>  +m  >  + 

S  <  ElLo^iW</’*(OEr=o^fcW<^fc(0  Ez=o^;W0KO,<(>j  >  + 

^  <  mi=oXi{t)^i{0  T.l=oXk{t)4>k{i)  TJi=QXi{t)(t)i{i),4>j  >=o 

/or  j=0,...,P. 


This  equation  can  be  written  as 


E*=o(^'*(^)  +  ta;i(t))  <  (j)i,(j)j  >  + 

^EjloEr=oE?Lo^*W^fcW^K^)  <  >  + 

^  Eto  Er=o  (Pj  >=  0 

/or  j=0,...,P. 


(29) 


This  system  of  j  equations  along  with  the  initial  conditions  needs  to  be  solved 
to  find  Xk{t)  in  Equation  (25).  Once  this  is  determined,  the  transformation  can  be 
approximated  the  same  way  as  in  the  linear  case. 

From  the  parameters  given  for  the  model  ODE  in  Table  (4),  we  determine 
the  approximate  transformation  for  the  nonlinear  problem.  For  this  example,  the 
stochastic  projection  method  described  uses  the  one- dimensional  Hermite  polynomials 
in  terms  of  normally  distributed  random  variables  (^)  with  a  mean  of  0  and  variance 
of  1  as  described  in  [4]. 

For  the  function  we  used  the  scaled  Hermite  polynomial 


(30) 


and  for  the  weight  function  w{^),  we  use  the  pdf  of  the  standard  normal  distribution 


28 


(31) 


The  limits  of  integration  used  for  the  inner  product  calculations  are  {—oo,  oo)  and 
i,j  =  0,...,P.  The  approximated  and  actual  transformations  are  shown  in  Fig¬ 
ure  (19)  where  g{^)  is  the  solid  plot  and  g{^)  is  the  dashed  plot.  The  Mathematica 
code  for  the  Hermite-chaos  method  can  be  found  in  Appendix  B.l. 


Figure  19:  Nonlinear  ODE  Transformation  Using  Hermite-chaos  of  Degree  2  and  4. 


3.2.2  Approximating  the  Transformation  g-.X^Y  Using  the  Uni¬ 
form  Sampling  Method.  The  uniform  sampling  method  used  for  the  nonlinear 
case  is  exactly  the  same  as  the  linear  case.  We  approximate  the  transformation  g{x) 
using  a  set  {{xi,g{xi)}  where  {xi, . . .  ,Xn}  is  not  random  but  is  chosen  to  be  evenly 
spaced  on  the  interval  [a,  h].  Then,  {g{xi), . . . ,  g{xn)}  is  calculated  by  solving  the  ode 
system  for  each  Xi.  Using  this  data,  an  approximation  g{x)  ~  g{x)  is  constructed.  The 
approximated  transformation  is  shown  in  Figure  (20)  where  the  actual  transformation 
is  the  solid  plot  and  the  approximated  transformation  is  the  dashed  plot.  These  two 
plots  are  almost  identical.  The  Mathematica  code  for  the  uniform  sampling  method 
can  be  found  in  Appendix  B.2.  Once  we  have  g{X),  the  transformation  method  can 
be  used  to  obtain  the  cdf.  Then  by  differentiating  the  cdf,  the  pdf,  in  Figure  (21), 
is  found.  This  pdf  obtained  with  the  uniform  sampling  method  can  be  compared  to 
the  pdf  found  using  Monte  Carlo  Simulation.  The  Monte  Carlo  pdf  is  also  shown  in 
Figure  (21). 
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Figure  20:  Nonlinear  ODE  Transformation  Using  Uniform  Sampling  Method  With 
256  Samples 


Figure  21:  PDF  of  Nonlinear  ODE  Using  Uniform  Sampling  Method  With  256 
Samples  and  Monte  Carlo  Simulation  With  50,000  Runs 
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IV.  Comparison  of  Hermite-chaos  and  Uniform  Sampling 

Method  Used  to  Approximate  the  Transformation 

As  demonstrated  previously,  Hermite-chaos  and  the  uniform  sampling  method  both 
worked  well  to  approximate  the  transformation  for  linear  and  nonlinear  ODE’s  when 
the  functional  form  of  g{X)  was  not  explicitly  known.  This  chapter  will  provide  a 
comparison  of  these  methods  based  on  the  accuracy  of  the  approximation  at  different 
levels  and  the  amount  of  computational  cost  required  for  each  method. 

4-1  Accuracy  and  Computational  Cost  Evaluation  of  Approximated  Trans¬ 
formation 

The  determination  of  how  accurate  the  approximations  are  will  be  made  by 
computing  the  weighted  error  for  different  levels  of  each  approximation  method.  The 
method  used  to  assess  the  accuracy  is 


\T{i)-f{i)yw{Odx  (32) 

-6 

where  T(^)  is  the  actual  transformation,  T(^)  is  the  approximate  transformation, 
and  w(0  is  the  pdf  of  the  standard  normal,  iV(0, 1),  distribution  since  we  assume 
the  input  random  variable  can  be  written  in  terms  of  the  standard  normal  random 
variable.  The  Mathematica  code  for  the  weighted  error  calculations  can  be  found  in 
Appendix  A.  The  computational  cost  was  measured  by  the  time,  in  seconds,  each 
method  took  to  execute  using  Mathematica.  Only  the  amount  of  computer  code 
necessary  to  obtain  the  approximated  transformation  was  processed.  The  TimeUsed 
function  in  Mathematica  recorded  the  time  spent  in  the  kernel  to  execute  the  selected 
code.  This  is  the  method  for  recording  the  computational  time  that  will  be  used 
throughout  this  paper. 
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4-1.1  Linear  ODE.  The  linear  ODE  under  analysis  is  the  model  problem 
(1)  with  the  parameters  in  Table  (4)  from  Chapter  3. 

4.1.1.1  Hermite- chaos.  To  evaluate  the  error,  when  using  Hermite- 
chaos  to  hnd  the  approximation,  the  weighted  error  was  calculated  using  Hermite 
polynomials  with  degree  2,  4,  6,  8,  10,  and  12.  As  the  degree  of  the  Hermite  polyno¬ 
mial  increases,  the  accuracy  of  g{X)  increases  dramatically.  However,  the  computa¬ 
tional  cost  also  increases  signihcantly.  The  graphs  in  Figure  (22)  show  the  increase 
in  accuracy  between  hermite  polynomial  of  degree  2  and  4  respectively. 

g  g 


Figure  22:  Approximate  (dashed)  and  Actual  (solid)  Transformation  oiY  =  g{X) 
Using  Hermite-chaos  of  Degree  2  and  4  for  Linear  ODE. 

The  weighted  error  for  degree  2  and  4  are  1.10399  and  0.306322  respectively. 
The  error  continues  to  decrease  for  each  increase  in  even  degree.  The  log-log  graph 
in  Figure  (23)  shows  the  weighted  error  for  each  degree  of  the  hermite  polynomial 
analyzed.  We  see  that  the  error  decrease  is  exponential. 

4- 1.1. 2  Uniform  Sampling  Method.  For  this  method,  the  weighted 
error  was  calculated  using  increasing  number  of  samples  (8,  16,  32,  64,  128,  256,  512, 
1024,  2048,  and  4096).  These  are  simply  determined  by  letting  the  sample  size  be  2” 
where  n  =  3, 4,  5  ...  12.  As  the  number  of  samples  increases,  the  accuracy  of  g{X)  also 
increases.  The  graphs  in  Figure  (24)  show  the  increase  in  accuracy  between  samples 
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Weighted  Error 


Figure  23:  Log-Log  Weighted  Error  Plot  by  Degree  of  Hermite  Polynomial  for  Lin¬ 
ear  ODE. 

of  32  and  64  respectively.  The  weighted  error  for  32  and  64  samples  are  0.0389027  and 
0.00973277  respectively.  The  error  continues  to  decrease  for  each  increase  in  number 
of  samples.  The  log-log  graph  in  Figure  (25)  shows  the  weighted  error  for  each  sample 
set  analyzed. 


Figure  24:  Approximate  (dashed)  and  Actual  (solid)  Transformation  oiY  =  g{X) 
Using  Uniform  Sampling  Method  With  32  and  64  Samples  for  Linear  ODE. 


In  Table  (5),  the  accuracy  and  computational  time  of  the  uniform  sampling 
method  and  Hermite-chaos  are  compared.  We  see  that  the  uniform  sampling  method 
produces  better  results  more  efficiently.  When  comparing  the  Hermite-chaos  method 
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Weighted  Error 


Figure  25:  Log-Log  Weighted  Error  Plot  by  Number  of  Samples  for  Linear  ODE. 


of  degree  4  and  the  uniform  sampling  method  with  64  samples,  we  see  that  the  uniform 
sampling  method  produces  a  more  accurate  result  in  roughly  the  same  amount  of 
computational  time.  This  is  also  the  case  for  Hermite-chaos  of  degree  8  and  256 
samples.  This  difference  disappears  as  polynomial  degree  and  number  of  samples 
increase. 


Hermite-Chaos 

Uniform  Sampling 

Degree 

Time(Seconds) 

Error 

Samples 

Time(Seconds) 

Error 

2 

0.733 

1.10399 

32 

.687 

0.0389027 

4 

0.876 

0.306322 

64 

.860 

0.00973277 

6 

1.219 

0.104656 

128 

1.11 

0.00243358 

8 

1.656 

0.007349 

256 

1.671 

0.000601922 

10 

2.376 

0.0001943 

512 

2.751 

0.000150971 

Table  5:  Comparison  of  Computational  Time  and  Error  Between  Hermite- 
Chaos  and  Uniform  Sampling  Method  Used  to  Approximate  the  Transfor¬ 
mation  for  the  Linear  ODE 
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When  analyzing  the  error  rate  of  convergence,  for  the  nniform  sampling  method, 
we  see  that  the  error  converges  at  a  qnadratic  rate  which  is  consistent  with  linear 
interpolation  methods.  The  error  rate  of  convergence  (ROC)  for  the  linear  case  was 
calcnlated  nsing  the  following  formnla 


ROC 


]  ( Errori20^\ 

y  V  Error{4^096)  ^ 


log{ 


4096 

2048 


) 


(33) 


The  weighted  errors  for  2048  and  4096  samples  are  8.8297x10  ®  and  2.2466x10  ® 
respectively.  The  error  rate  of  convergence  was  calculated  to  be  1.97464. 

4.1.2  Nonlinear  ODE.  The  procedures  for  assessing  the  weighted  error 
of  the  nonlinear  ODE  are  exactly  the  same  as  the  linear  ODE.  The  ’’actual  transfor¬ 
mation”  was  obtained  by  the  sampling  method  with  8192  samples.  This  will  be  used 
to  compare  the  accuracy  of  the  different  levels  of  each  approximation  method.  The 
Mathematica  code  for  the  weighted  error  calculations  can  be  found  in  Appendix  B. 
The  following  results  were  obtained. 

4. 1.2.1  Hermite- chaos.  Again,  as  the  degree  of  the  hermite  poly¬ 
nomial  increased,  the  accuracy  of  g{X)  increased  dramatically.  The  graphs  in  Fig¬ 
ure  (26)  show  the  increase  in  accuracy  between  hermite  polynomial  of  degree  2  and 
4  respectively.  The  weighted  error  for  degree  2  and  4  are  0.167758  and  0.0164006 
respectively.  The  error  continues  to  decrease  for  each  increase  in  even  degree.  The 
log-log  graph  in  Figure  (27)  shows  the  weighted  error  for  each  degree  of  the  hermite 
polynomial  analyzed.  Again,  we  see  that  the  error  decrease  is  exponential. 
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Figure  26:  Approximate  (dashed)  and  Actual  (solid)  Transformation  oiY  =  g{X) 
Using  Hermite-chaos  of  Degree  2  and  4  for  Nonlinear  ODE. 


Weighted  Error 


Figure  27:  Log-Log  Weighted  Error  Plot  by  Degree  of  Hermite  Polynomial  for  Non¬ 
linear  ODE. 

4- 1.2. 2  Uniform  Sampling  Method.  In  the  nonlinear  case,  just  as 
the  linear,  when  the  number  of  samples  increases  the  accuracy  of  g{X)  increases.  The 
graphs  in  Figure  (28)  show  the  increase  in  accuracy  between  samples  of  32  and  64. 
The  weighted  error  for  32  and  64  samples  are  0.00802643  and  0.00200724  respectively. 
The  error  continues  to  decrease  for  each  increase  in  number  of  samples.  The  graph  in 
Figure  (29)  shows  the  weighted  error  for  each  sample  set  analyzed. 

In  Table  (6),  the  accuracy  and  computational  time  of  the  nniform  sampling 
method  and  Hermite-chaos  are  compared  for  this  nonlinear  example.  We  see  that 
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Figure  28:  Approximate  (dashed)  and  Actual  (solid)  Transformation  oiY  =  g{X) 
Using  Uniform  Sampling  Method  With  32  and  64  Samples  for  Nonlinear  ODE. 


Weighted  Error 


Figure  29:  Log-Log  Weighted  Error  Plot  by  Number  of  Samples  for  Nonlinear  ODE. 

the  uniform  sampling  method  produces  better  results  more  efficiently.  When  com¬ 
paring  the  Hermite-chaos  method  of  degree  4  and  the  uniform  sampling  method  with 
32  samples,  we  see  that  the  uniform  sampling  method  produces  a  more  accurate  re¬ 
sult  in  much  less  computational  time.  There  is  a  more  dramatic  difference  between 
the  Hermite-chaos  of  degree  8  and  256  samples.  While  the  error  is  approximately 
the  same,  the  Hermite-chaos  method  to  approximate  the  transformation  takes  close 
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to  840  times  longer  than  the  sampling  method.  In  contrast  to  the  linear  example 
described  previonsly,  the  nniform  sampling  method  ontperforms  the  Hermite-chaos 
method  significantly  for  the  nonlinear  example.  The  compntational  time  reqnired  to 
rnn  the  Hermite-chaos  method  is  mostly  spent  on  the  calcnlation  of  the  inner  prodncts 
in  Eqnation  (29).  When  these  inner  prodncts  were  precalcnlated  and  the  resnlts  were 
snbstitnted  into  the  Mathematica  code,  the  compntational  cost  was  greatly  rednced. 
Once  the  needed  inner  prodncts  are  calculated,  they  can  be  substituted  in  for  any  de¬ 
gree  of  Hermite  polynomial  needed  for  finding  g{X)  from  the  model  problem  in  (24). 
The  limitation  on  this  procedure  is  that  the  inner  products  have  to  be  recalculated  for 
different  types  of  nonlinearities.  The  results  using  the  precalculated  inner  products 
are  shown  in  Table  (7). 


Hermite-Chaos 

Sampling  Method 

Degree 

Time(Seconds) 

Error 

Samples 

Time(Seconds) 

Error 

2 

12.296 

0.1677582 

32 

0.780 

0.00802643 

4 

136.86 

0.0164006 

64 

1.11 

0.00200724 

6 

632.20 

0.0013722 

128 

1.609 

0.000501754 

8 

2265.7 

0.0001269 

256 

2.688 

0.000125043 

10 

7114.8 

0.0000145 

512 

4.891 

0.0000275471 

Table  6:  Comparison  of  Computational  Time  and  Error  Between  Hermite- 
Chaos  and  Uniform  Sampling  Method  Used  to  Approximate  the  Transfor¬ 
mation  for  the  Nonlinear  ODE 


Hermite-Chaos 

Degree 

Time(Seconds) 

Error 

2 

1.736 

0.16776 

4 

1.797 

0.01640 

6 

1.969 

0.001372 

8 

2.282 

0.000127 

10 

3 

0.0000145 

Table  7:  Computational  Time  and  Error  of  Hermite-Chaos  Method  With 
Precalculated  Inner  Products  for  the  Nonlinear  ODE 
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When  analyzing  the  error  rate  of  convergence,  for  the  nniform  sampling  method, 
we  see  that  the  error  converges  at  a  qnadratic  rate  which  is  consistent  with  linear 
interpolation  methods.  The  error  rate  of  convergence  for  the  nonlinear  case  was 
calcnlated  nsing  the  following  formula 


ROC 


]  ( Errori20^\ 

y  V  Error{4^096)  ^ 


log{ 


4096 

2048 


) 


(34) 


The  weighted  errors  for  2048  and  4096  samples  are  1.6645x10  ®  and  3.9598x10  ^ 
respectively.  The  error  rate  of  convergence  was  calculated  to  be  2.07156. 

4.2  Accuracy  Evaluation  of  Approximated  PDF 

Now  that  the  accuracy  of  g,  approximated  with  Hermite-chaos  and  the  uniform 
sampling  method,  has  been  analyzed,  we  will  look  at  the  precision  of  the  pdf  generated 
from  these  approximated  transformations.  As  discussed  in  Chapter  2,  there  were 
two  methods  of  approximating  the  pdf;  the  transformation  method  and  Monte  Carlo 
simulation.  The  approximated  pdfs  will  be  compared  to  the  pdfs  obtained  from  the 
actual  transformations  for  the  linear  and  nonlinear  examples.  The  method  used  to 
assess  the  accuracy  is 


-?<*/(!/)) Vj/ 


(35) 


where  pdf{y)  is  the  actual  pdf  and  pdf{y)  is  the  approximated  pdf.  The  interval 
for  the  linear  example  will  be  (a,  6)  =  (—3,4.1)  while  the  interval  for  the  nonlinear 
example  will  be  (a,  6)  =  (-10,-4).  These  modified  intervals  are  necessary  due  to 
erratic  behavior  at  the  edges  of  the  pdf’s  from  the  use  of  the  interpolation  function 
in  Mathematica.  The  limits  selected  adequately  capture  where  the  actual  pdf  has 
positive  probability. 
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4.2.1  Linear  ODE. 


4.2. 1.1  Transformation  Method.  To  evaluate  the  error  when  find¬ 
ing  the  pdf  with  the  transformation  method  using  Hermite-chaos,  the  error  was  cal¬ 
culated  using  the  actual  pdf  and  the  approximated  pdf  using  hermite  polynomials 
with  degree  2,  4,  6,  8,  10,  and  12.  The  graphs  in  Figure  (30)  show  the  increase  in 
accuracy  between  hermite  polynomials  of  degree  2  and  4.  The  error  for  degree  2  and 
4  are  0.0228001  and  0.0154914  respectively.  The  error  continues  to  decrease  for  each 
increase  in  even  degree.  The  graph  in  Figure  (31)  shows  the  error  for  each  degree  of 
the  hermite  polynomial  analyzed. 


fy 
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Figure  30:  Approximate  (solid)  and  Actual  (dashed)  pdf{X)  (Via  Transformation 
Method)  Using  Hermite-chaos  of  Degree  2  and  4  for  the  Linear  ODE. 


To  evaluate  the  error  when  finding  the  pdf  with  the  transformation  method  using 
uniform  sampling,  the  error  was  calculated  using  the  actual  pdf  and  the  approximate 
pdf  with  sample  sizes  of  2”  where  n  =  4,  5, ...  12.  The  sample  size  of  8  was  too 
small  to  compute  a  pdf  via  the  transformation  method.  The  graphs  in  Figure  (32) 
show  the  increase  in  accuracy  between  samples  of  32  and  64.  The  error  for  32  and 
64  samples  are  0.0640926  and  0.0178453  respectively.  The  error  continues  to  decrease 
for  each  increase  in  number  of  samples.  The  graph  in  Figure  (33)  shows  the  error  for 
each  sample  set  analyzed.  The  comparison  between  Hermite-chaos  and  the  uniform 
sampling  method,  in  terms  of  pdf  error  is  shown  in  Table  (8). 
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Figure  31:  Log-Log  Error  Plot  for  pdf{X)  (Via  Transformation  Method)  by  Degree 
of  Hermite  Polynomial  for  the  Linear  ODE. 


Figure  32:  Approximate  (solid)  and  Actual  (dashed)  pdf{X)  (Via  Transformation 
Method)  Using  Sampling  Method  With  32  and  64  Samples  for  the  Linear  ODE. 


4 -2. 1.2  Monte  Carlo  Method.  To  obtain  a  pdf  from  the  histogram 
generated  by  Monte  Carlo  simulation,  the  midpoints  of  the  frequency  rectangles  were 
linearly  interpolated  to  form  the  pdf  function  to  be  analyzed.  When  using  Hermite- 
chaos  to  approximate  the  transformation,  the  error  was  calculated  using  the  actual 
pdf  and  the  approximated  pdf  found  using  hermite  polynomials  with  degree  2,  4, 
6,  8,  10,  and  12.  The  graphs  in  Figure  (34)  show  the  increase  in  accuracy  between 
hermite  polynomial  of  degree  2  and  4.  The  error  for  degree  2  and  4  are  0.423889 
and  0.0323146  respectively.  The  error  continues  to  decrease  for  each  increase  in  even 
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Figure  33:  Log-Log  Error  Plot  for  pdf{X)  (Via  Transformation  Method)  by  Number 
of  Samples  for  the  Linear  ODE. 


Hermite-Chaos 

Sampling  Method 

Degree 

Error 

Samples 

Error 

2 

285.081 

32 

0.0640926 

4 

0.0228001 

64 

0.0178453 

6 

0.0154914 

128 

0.00886257 

8 

0.00207821 

256 

0.00595147 

10 

0.000181367 

512 

0.00254538 

12 

.000003249 

4096 

.000001339 

Table  8:  Comparison  of  Error  Between  Hermite-chaos  and  Uniform 
Sampling  Method  Used  to  Approximate  the  PDF  (Via  Transformation 
Method)  of  the  Linear  ODE 


degree.  The  graph  in  Figure  (35)  shows  the  error  for  each  degree  of  the  hermite 
polynomial  analyzed. 

To  evaluate  the  error  when  finding  the  pdf  with  Monte  Carlo  simulation  using 
the  uniform  sampling  method,  the  error  was  calculated  using  the  actual  pdf  and 
the  approximate  pdf  with  sample  sizes  of  2”  where  n  =  3,5,. ..12.  The  graphs  in 
Figure  (36)  show  the  increase  in  accuracy  between  samples  of  32  and  64  respectively. 
The  error  for  32  and  64  samples  are  0.0885416  and  0.024061  respectively.  The  error 
fluctuates  between  each  number  of  samples.  This  is  typical  of  Monte  Carlo  simulation 
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Figure  34:  Approximate  (solid)  and  Actual  (dashed)  pdf{X)  (Via  Monte  Carlo 
Method)  Using  Hermite-chaos  of  Degree  2  and  4  for  the  Linear  ODE. 


Error 


Figure  35:  Log-Log  Error  Plot  for  pdf{X)  (Via  Monte  Carlo  Method)  by  Degree  of 
Hermite  Polynomial  for  the  Linear  ODE. 


due  to  the  random  input  used  to  obtain  the  pdf.  The  graph  in  Figure  (37)  shows  the 
error  for  each  sample  set  analyzed. 
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Figure  36:  Approximate  (solid)  and  Actual  (dashed)  pdf{X)  (Via  Monte  Carlo 
Method)Using  Uniform  Sampling  With  32  and  64  Samples  for  the  Linear  ODE. 


Error 


Figure  37:  Log  Log  Error  Plot  for  pdf{X)  (Via  Monte  Carlo  Method)  by  Number 
of  Samples  for  the  Linear  ODE. 

4-2.2  Nonlinear  ODE. 

4-2.2. 1  Transformation  Method.  When  using  Hermite-chaos,  the 
error  was  calculated  using  the  actual  pdf  and  the  approximated  pdf  using  hermite 
polynomials  with  degree  2,  4,  6,  8,  10,  and  12.  This  is  the  same  procedure  used  for 
the  linear  example.  The  graphs  in  Figure  (42)  show  the  increase  in  accuracy  between 
hermite  polynomial  of  degree  2  and  4  respectively.  The  error  for  degree  2  and  4  are 
71.199  and  1.09511  respectively.  The  error  continues  to  decrease  for  each  increase  in 
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even  degree  until  after  degree  6  where  the  error  bottoms  out.  The  graph  in  Figure  (39) 
shows  the  error  for  each  degree  of  the  hermite  polynomial  analyzed. 

fy  fy 


Figure  38:  Approximate  (Solid)  and  Actual  (dashed)  pdf{X)  (Via  Transformation 
Method)  Using  Hermite-chaos  of  Degree  2  and  4  for  Nonlinear  ODE. 


Error 


Figure  39:  Error  Plot  for  pdf{X)  (Via  Transformation  Method)  by  Degree  of  Her¬ 
mite  Polynomial  for  Nonlinear  ODE. 

When  using  uniform  sampling,  the  error  was  calculated  using  the  actual  pdf  and 
the  approximated  pdf  with  sample  sizes  of  2"^  where  n  =  3,5,..  .12.  The  graphs  in 
Figure  (40)  show  the  increase  in  accuracy  between  samples  of  32  and  64  respectively. 
The  error  for  32  and  64  samples  are  0.723152  and  0.705609  respectively.  The  error 
continues  to  decrease  for  each  increase  in  number  of  samples.  The  graph  in  Figure  (41) 
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shows  the  weighted  error  for  each  sample  set  analyzed.  The  comparison  between 
Hermite-chaos  and  the  sampling  method,  in  terms  of  pdf  error  is  shown  in  Table  (9). 
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Fignre  40:  Approximate  (solid)  and  Actnal  (dashed)  pdf{X)  (Via  Transformation 
Method)  Using  Uniform  Sampling  Method  With  32  and  64  Samples  for  the  Nonlinear 


ODE. 


Error 


Figure  41:  Error  Plot  for  pdf{X)  (Via  Transformation  Method)  by  Number  of  Sam¬ 
ples  for  Nonlinear  ODE. 


4 -2. 2. 2  Monte  Carlo  Method.  When  using  Hermite-chaos,  the 
error  was  calculated  using  the  actual  pdf  and  the  approximated  pdf  using  hermite 
polynomials  with  degree  2,  4,  6,  8,  10,  and  12.  This  is  the  same  procedure  used  for 
the  linear  example.  The  graphs  in  Figure  (42)  show  the  increase  in  accuracy  between 
hermite  polynomial  of  degree  2  and  4.  The  error  for  degree  2  and  4  are  0.903489  and 


46 


Hermite-Chaos 

Sampling  Method 

Degree 

Error 

Samples 

Error 

2 

71.199 

32 

0.723152 

4 

1.09511 

64 

0.705609 

6 

0.655017 

128 

0.68144 

8 

0.615786 

256 

0.612364 

10 

0.612495 

512 

0.589525 

12 

0.612401 

4096 

0.380311 

Table  9:  Comparison  of  Error  Between  Hermite-chaos  and  Uniform 
Sampling  Method  Used  to  Approximate  the  PDF  (Via  Transformation 
Method)  of  the  Nonlinear  ODE 


0.758707  respectively.  The  error  decreases  from  degree  2  to  4  and  fluctuates  around 
0.758  for  degree  greater  than  4.  The  graph  in  Figure  (43)  shows  the  error  for  each 
degree  of  the  hermite  polynomial  analyzed. 


Figure  42:  Approximate  (solid)  and  Actual  (dashed)  pdf{X)  (Via  Monte  Carlo  Sim¬ 
ulation)  Using  Hermite-chaos  of  Degree  2  and  4  for  the  Nonlinear  ODE. 


To  evaluate  the  error  when  finding  the  pdf  with  Monte  Carlo  simulation  using 
the  uniform  sampling  method,  the  error  was  calculated  using  the  actual  pdf  and 
the  approximated  pdf  with  sample  sizes  of  2"^  where  n  =  3,5,..  .12.  The  graphs  in 
Figure  (44),  where  pdf{X)  is  the  solid  plot  and  pdf{X)  is  the  dashed  plot,  show  the 
increase  in  accuracy  between  samples  of  32  and  64.  The  errors  for  32  and  64  samples 
are  0.761875  and  0.758375  respectively.  The  error  decreases  for  each  increase  in  level 
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Figure  43:  Error  Plot  for  pdf{X)  (Via  Monte  Carlo  Simulation)  by  Degree  of  Her- 
mite  Polynomial  for  the  Nonlinear  ODE. 


of  samples  up  to  128.  From  there  the  error  fluctuates  above  and  below  0.757.  The 
graph  in  Figure  (45)  shows  the  error  for  each  sample  set  analyzed. 
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Figure  44:  Approximate  (solid)  and  Actual  (dashed)  pdf{X)  (Via  Monte  Carlo  Sim¬ 
ulation)  Using  Uniform  Sampling  Method  With  32  and  64  Samples  for  the  Nonlinear 
ODE. 


It  is  important  to  point  out  that  the  pdfs  generated  by  Monte  Carlo  simulation 
have  many  oscillations  making  it  difficult  to  reduce  the  error.  There  are  methods 
available  that  smooth  out  the  pdf,  but  these  are  extra  operations  that  need  to  be  per¬ 
formed.  When  using  the  transformation  method,  there  are  no  extra  steps  involved  to 
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Figure  45:  Error  Plot  for  pdf{X)  (Via  Monte  Carlo  Simulation)  by  Number  of 
Samples  for  the  Nonlinear  ODE. 


get  a  meaningful  pdf  function.  The  transformation  method  provides  a  more  accurate 
representation  of  the  pdf  generated  from  the  actual  or  approximated  transformation. 
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V.  Nonuniform  Sampling  Method 

In  the  examples  analyzed  in  the  last  chapter,  uniform  sampling  was  shown  to  be 
more  accurate  and  faster  than  the  traditional  Monte  Carlo  and  stochastic  projection 
methods  currently  used  in  uncertainty  analysis.  We  extend  the  idea  further  and 
consider  a  nonuniform  sampling  approach.  This  idea  stems  from  the  notion  that 
the  input  random  variable  is  from  a  standard  normal  distribution.  In  this  case,  the 
values  of  the  random  input  that  are  closer  to  the  mean  have  a  higher  probability  of 
occurrence  than  the  values  that  are  2  or  3  standard  deviations  away  from  the  mean. 
This  implies  that  a  collection  of  samples  that  are  a  more  realistic  representation  of 
the  actual  input  values  would  produce  a  more  accurate  transformation  and,  therefore, 
produce  a  better  pdf  of  the  output  random  variable. 

To  produce  a  nonuniform  sample  based  on  the  A^(0, 1)  distribution  of  the  input 
random  variable,  we  can  generate  a  set  of  evenly  spaced  nodes  on  [—a,  a],  (denoted 
as  N),  and  then  evaluate  the  standard  normal  CDF,  shown  in  Figure  (46),  at  these 
nodes  (denoted  as  C).  Then  the  pairs,  {(cj,nj)  :  c*  e  C, n*  e  N,for  i  =  l,...,m}, 
are  interpolated,  producing  a  graph  of  the  interpolated  inverse  CDF.  We  evaluate  the 
inverse  CDF  to  produce  a  sample  from  the  Normal[0,l]  distribution.  To  illustrate 
this  method,  we  generate  the  inverse  CDF,  shown  in  Figure  (46),  using  500  equally 
spaced  nodes  on  [—6,6].  We  then  evaluate  this  inverse  CDF  at  at  each  node  of  N 
where  N  =  {1/16, 1/8,  3/16, 1/4, ...,  13/16,  7/8, 15/16}.  The  sample  generated  is 
shown  in  Figure  (47).  In  addition  to  these  samples,  the  endpoints  of  the  interval 
[—a,  a]  and  two  additional  points  are  added  to  the  set  in  order  to  capture  those  areas 
of  the  transformation. 

The  nonuniform  sampling  method  was  tested  on  the  nonlinear  ODE  model  prob¬ 
lem  (24)  from  Chapter  3.  Once  the  nonuniform  sample  was  generated,  the  same  algo¬ 
rithm  to  obtain  the  approximated  transformation  was  used  as  in  the  uniform  sampling 
method.  Just  as  the  uniform  sampling  case,  when  the  number  of  samples  increase  the 
accuracy  of  g{X)  increases.  The  graphs  in  Figure  (48)  show  the  increase  in  accuracy 
between  samples  of  35  and  67.  The  weighted  error  for  35  and  67  samples  are  0.040645 


50 


Figure  46:  CDF  and  Inverse  CDF  of  N(0,1)  Distribution. 
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Figure  47:  Nonuniform  Sample  Generated  Based  on  N(0,1)  Distribution. 


and  0.0230532  respectively.  The  graph  in  Figure  (49)  shows  the  weighted  error  for 
each  sample  set  analyzed.  The  Mathematica  code  used  for  the  nonuniform  sampling 
method  can  be  found  in  Appendix  B.3. 


Figure  48:  Approximate  (dashed)  and  Actual  (solid)  Transformation  oiY  =  g{X) 
Using  Nonuniform  Sampling  Method  With  35  and  67  Samples  for  Nonlinear  ODE. 
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Weighted  Error 


Figure  49:  Log-Log  Weighted  Error  Plot  by  Number  of  Nonuniform  Samples  for 
Nonlinear  ODE. 

In  Table  (10),  the  accuracy  of  the  uniform  and  nonuniform  sampling  methods 
are  compared  for  this  nonlinear  example.  We  see  that  the  uniform  sampling  method 
produces  more  accurate  results  at  each  level  of  samples.  This  loss  of  accuracy  can 
be  attributed  to  the  lack  of  points  near  the  ends  of  the  interval  [—a,  a]  and  where 
the  second  derivative  of  the  transformation  is  non  zero.  These  improvements  could 
greatly  reduce  the  error,  but  will  not  be  covered  here.  Further  research  needs  to  be 
conducted  to  explore  and  improve  this  idea  of  nonuniform  sampling. 


Nonuniform  Sampling 

Uniform  Sampling 

Samples 

Error 

Samples 

Error 

35 

0.040645 

32 

0.00802643 

67 

0.0230532 

64 

0.00200724 

131 

0.0131352 

128 

0.000501754 

259 

0.0075023 

256 

0.000125043 

515 

0.0042869 

512 

0.0000275471 

Table  10:  Comparison  of  Weighted  Error  Between  Nonuniform  and  Uni¬ 
form  Sampling  Methods  Used  to  Approximate  the  Transformation  of  the 
Nonlinear  ODE 
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VI.  Conclusion  and  Further  Research 


Uncertainty  analysis  plays  an  important  role  in  engineering  design  and  analysis.  For  a 
long  time,  Monte  Carlo  simulations  and  stochastic  projection  methods  have  been  used 
to  determine  output  distributions  of  these  systems.  These  methods,  while  effective, 
can  be  cumbersome  and  inefficient. 

The  introduction  of  sampling  combined  with  the  transformation  method  can 
open  up  new  doors  in  determining  crucial  information  about  random  systems.  As 
shown  in  this  research,  the  uniform  sampling  method  was  an  effective  way  to  approx¬ 
imate  the  transformation  function  of  the  linear  and  nonlinear  ordinary  differential 
equation  model  problems.  This  method  is  also  quite  simple  to  execute.  In  chapter  3, 
the  methodology  for  the  Hermite-chaos  was  shown  for  the  linear  and  nonlinear  ODEs. 
The  calculations  required  were  complicated  and  consumed  a  great  amount  of  time  to 
complete.  The  uniform  sampling  method,  however,  was  simple  and  easy  to  calculate 
without  sacrificing  accuracy.  Monte  Carlo  simulation  was  also  simple  to  conduct,  but 
the  limited  accuracy  and  the  necessity  of  running  large  amounts  of  random  samples 
through  the  system  make  it  less  appealing.  By  using  the  transformation  method, 
finding  the  distribution  was  faster  and  more  accurate  than  Monte  Carlo  Simulation. 

While  the  uniform  sampling  and  transformation  method  outperformed  Hermite- 
chaos  and  Monte  Carlo  simulation  for  the  two  specific  examples  discussed  in  chapter 
3,  more  research  needs  to  be  conducted.  This  idea  can  be  extended  to  systems  with 
more  than  one  random  input.  The  uniform  sampling  method  should  also  be  tested 
on  real  world  applications.  The  notion  of  nonuniform  sampling  can  also  be  extended 
and  improved  to  validate  its  effectiveness.  In  addition,  finite  elements  could  be  used 
to  approximate  the  transformation  and  the  results  could  be  compared  to  sampling 
results. 
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The  sampling  and  transformation  method,  evaluated  in  this  research,  can  serve 
as  a  benchmark  for  future  improvements  in  the  held  of  uncertainty  analysis.  As  men¬ 
tioned  before,  many  advancements  can  be  made  on  this  very  simple  idea.  System 
reliability,  interaction  analysis,  characteristic  determination,  and  many  other  impor¬ 
tant  applications  can  be  simplihed  if  the  sampling  and  transformation  method  idea 
was  further  proven  to  be  an  effective  analytic  tool. 
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Appendix  A.  Mathematica  Code  for  the  Linear  ODE  to 
Approximate  the  Transformation 


A.l  Hermite- chaos  of  Degree  4 


<  <  Statistics'  ContinuousDistributions' 

<  <  Graphics'  Graphics' 

npdf[t_]  =  PDF[NormalDistribution[0, 

/[t_]  =  npdf[t]; 

Define  the  parameters  and  conditions  of  the  linear  ODE. 


a  =  —6; 
6  =  6; 
pd  =  4; 
m  =  10; 
kb  =  10; 
kt  =  2; 
c  =  4; 


V  =  25; 


HermiteH[j, 

y/¥if. 


Create  A  and  B  matrices  for  S. 


ipa  =  Identity  Matrix[pd  + 1]; 

MatrixFor  m  [ipa] ; 

ipb  =  Table[0,  {i,  0,  pd},  {j,  0,  pd}]; 

For[«  =  0, «  <  pd, «++, 

For\j  =  ij  <  pd,  j++, 

ipb[[i  +  l,j  +  1]]  =  Chop[NIntegrate[tp[t,i]p[t,  j]/[t],  {t,  —Infinity,  Infinity}]]; 
ipb[[7  +  !,*  +  !]]=  ipb[[i  +  1,  j  +  1]]; 

] 
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] 

MatrixFor  m  [ipb] ; 

md  =  Dimensions[ipa][[l]]; 

aall  =  Table[0,  {i,  1,  md},  {j,  1,  md}]; 
aal2  =  IdentityMatrix[md]; 

aa21  =  — ((^)IdentityMatrix[md]  +  (^)Inverse[ipa].ipb); 
aa22  =  —  ^IdentityMatrbcfmdj; 

aa  =  Chop[Join[Transpose[Join[aall,  aal2]],Transpose[Join[aa21,aa22]]]] 
MatrixForm[aa] ; 

Solve  the  system  of  linear  ODEs. 

tf  =  10; 

ua  =  Table  [M[i][t],  {i,  1, 2md}]; 
upa  =  Table[M[i]'[t],  {f,  1, 2md}]; 

ics  =  Flatten[{Table[M[i][0]  ==  0,  {i,  1,  md}],  wfmd  +  1][0]  ==  t;, 
Table[M[z][0]  ==  0,  {*,  md  +  2, 2md}]}]; 
eqns  =  Flatten[Table[upa[[«]]  ==  (aa.ua)  [[«]],  {«,  1, 2md}]]; 
eqns  =  Flatten[{eqns,ics}]; 

MatrixForm  [eqns] ; 

soln  =  NDSolve[eqns,  Table[u[i],  {i,  1, 2md}],  {t,  0,  tf}]; 
ut  =  Table  [Evaluate  [Table  [M[i][f],  {*,  1,  md}]/.soln],  {t,  0,  tf}]; 

Table  [Evaluate  [Table  [M[i][t],  {z,md  +  1, 2md}]/.soln],  {f,  0,tf}]; 
cc  =  Flatten[ut[[Dimensions[ut][[l]]]]] 

{-0.537323,  -2.02877,  0.514555,  0.941688,  -0.584251} 

Calculate  the  approximated  transformation. 

as[t_]  =  Sum[cc[[i]]p[f,i  —  1],  {i,  l,md}]; 

TimeUsedO 

Plot[as[f],{t,  — 6,6},AxesLabel  — >  {"x",  "ghat(x)"}] 

Display [" ghat l.eps",%,  "EPS"] 
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nd  =  2«  +  1; 

id  =  Table[a  +  (i  -  1)^,  {*,  1,  nd}]; 
od  =  Table[0,  {i,  1,  nd}]; 

For[i  =  l,i  <  nd,i++, 
od[[i]]  =  as[id[[i]]] 

] 


This  is  the  actual  transformation  for  the  model  problem. 

qq[s-.t-]  =  -  4(^)(kb  +  skt)]]; 

9[s-]  =  sqrt[(^)2-44)(kb+akt)]^^P[-^]^^^^[fSqrt[(^)"  -  4(i)(kb  +  skt)]]; 

Plot[q'[s],  {s,  —6, 6},  AxesLabel  — >  {"x",  "g(x)"}] 

Display["g.eps",%,  "EPS"] 

Plot  the  actual  and  approximate  transformations  and  calculate  the  weighted  error. 

Plot[{q'[f],as[f]},{t,— 6, 6},  AxesLabel  — >  "g"},PlotStyle  — »  {GrayLevel[0],Dashing[{.03}]}] 

Export["gandgliat4.eps",  %,  "EPS"] 

WeightedError  =  Sqrt[NIntegrate[(g[t]  -as[t])V[t],{t,-6, 6}]] 
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g 


0.306322 

A. 2  Uniform  Sampling 

<  <  Statistics'  ContinuousDistributions' 

<  <  Graphics'  Graphics' 

Define  the  parameters  and  conditions  of  the  linear  ODE. 

npdf[t_]  =  PDF[NormalDistribution[0, 

/[t_]  =  npdf[f] 
a  =  —6; 

6  =  6; 
m  =  10; 
kb  =  10; 
kt  =  2; 
c  =  4; 
w  =  25; 
tf  =  10; 


This  is  the  actnal  transformation  for  the  model  problem. 
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qqls-t-]  =  gqrt[(^)2_4(\)(i,b+,i^)]Exp[-^]Sinh[|Sqrt[(^)"  -  4(i)(kb  +  skt)]] 

9W  =  sqrt[(^)^-44)(kb+.kt)]^^P[-^]S^^h[fSqrt[(^)"  -  4(^)(kb  +  skt)]] 

50e-*/5Sinh[|-y/^-|{10+2s)t] 

50Sinh[5^^-|(10+2s)] 

^W^5-hW+2s) 

Sample  the  function  defining  the  transformation  from  ^  to  x(T,^)  .  Use  this  set  of 
samples  to  approximate  the  transformation. 

Using  8  Samples, 

n  =  3; 

nnd  =  2”  +  1; 
h  =  _t2_- 

nd  =  Table[a  +  hk,  {k,  0,  nnd  —  1}]; 
odataS  =  Table[0,  {*,  1,  nnd}]; 

odataS  =  Flatten[Table[yl[a;]/.NDSolve[{yl'[t]  ==  y2[t], 

y2-|(]  ==  _(!±±«)yl|t]  _  ^y2|(|, 

yl|0]  ==  0.  y2|0|  ==  »}.  {yl.  y2},  {«,  0,  tt)y.x  tf,  {!.  1,  mid}]|; 

Interpolate  the  points  generated  from  the  uniform  sampling  method. 

T8  =  Ttanspose[{nd,  odataS}]; 

pp8  =  Interpolation[T8,InterpolationOrder  — »  1]; 

Continue  above  steps  for  each  level  of  samples.  Plot  the  actual  and  approximate 
transformations  and  calculate  the  weighted  error  for  each  level  of  samples. 

Using  8  Samples, 

Plot[{pp8[t],  g[t]},  {t,  -6, 6}] 

WES  =  Sqrt[NIntegrate[(g[t]  —  pp8[t])^/[t],  {t,  —6, 6}]] 
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0.601852 

Using  16  Samples, 

Plot[{ppl6[t], g[t]},  {t,  —6, 6},  AxesLabel  {"x",  "g"}, 
PlotStyle  ^  {Dashing[{.03}],  GrayLevel[0]}] 

WE16  =  Sqrt[NIntegrate[(q'[t]  —  ppl6[t])^/[t],  {t,  —6,6}]] 


g 


0.155178 

Repeat  this  process  for  the  rest  of  the  sample  levels.  Produce  the  error  plot  of  ap¬ 
proximated  transformation  for  the  linear  ODE. 

LogLogListPlot[{{8,  WE8},  {16,  WE16},  {32,  WE32},  {64,  WE64},  {128,  WE128},  {256,  WE256}, 
{512,  WE512},  {1024,  WE1024},  {2048,  WE2048},  {4096,  WE4096}}, 

AxesLabel  — >  {"Number  of  Samples",  "Weighted  Error"},  Plot  Joined  — >  True, 

Ticks  {{8, 32, 256, 1024, 4096},  {.005,  .025,  .15,  .0005,  .00005,  .000002}}] 
Export["linsamplingerror.eps",  %,  "EPS"] 


60 


Weighted  Error 


0 . 15 
0 . 025 
0 .005 

0 .0005 
0 . 00005 

2 . • 10“® 


of 


Samples 
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Appendix  B.  Mathematica  Code  for  the  Nonlinear  ODE  to 
Approximate  the  Transformation 

B.l  Hermite-chaos  of  Degree  4 

Define  the  parameters  and  conditions  of  the  linear  ODE. 

<  <  Statistics'  ContinuousDistributions' 

<  <  Graphics'  Graphics' 
a  =  —Infinity; 

b  =  Infinity; 
kb  =  1.; 
kt  =  0.02; 
xO  =  10.; 
tf  =  10; 
m  =  10.; 
kl  =  10.; 
k3b  =  1; 
k3t  =  0.02; 
w  =  0; 

npdf[t_]  =  PDF[NormalDistribution[0,  l],t]; 

/[t_]  =  npdf[t] 

.  ,  HermiteH[j,^] 

e~^ 

Calculate  the  inner  products. 

md  =  5; 

ipl  =  Chop[Table[iV[Integrate|j}[t,  i  —  l]p[t,j  —  l]p[t,  k  —  l]p[t,  I  —  l]/[t],  {t,  o,  6}],  30],  {«,  1,  md}, 
{j,  1,  md},  {k,  1,  md},  {Z,  1,  md}]]; 

ip2  =  Ghop[Table[iV[Integrate[tp[t,  i  —  l]p[t,j  —  l]p[t,  k  —  l\p[t,  I  —  1]/[Z],  {t,  a,  6}],  30],  {*,  1,  md}, 
if  1,  md},  {k,  1,  md},  {Z,  1,  md}]]; 
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Solve  the  system  of  nonlinear  ODEs. 

deqns  =  Chop  [Table  [M[s]”[t]  +  ^M[s][t]+ 

^Sum[M[i][t]M[7][t]M[fc][t]ipl[[i,  j,  k,  s]],  {i,  1,  md},  {j,  1,  md},  {k,  1,  md}]+ 
^Sum[M[i][t]M|j][t]M[fc][t]ip2[[i,  j,  k,  s]],  {i,  1,  md},  {j,  l,md},  {k,  1,  md}]  ==  0,  {s,  1,  md}]]; 
ics  =  Flatten[{M[l]  [0]  ==  xO,  Table[M[z][0]  ==  0.,  {i,  2,  md}],  Table[M[i]'[0]  ==  0.,  {i,  1,  md}]}]; 
eqns  =  Flatten[{deqns,ics}]; 

MatrbcForm  [eqns] ; 

soln  =  NDSolve[eqns,  Table[M[«],  {«,  1,  md}],  {t,  0,  tf}.  Method  ^  Adams]; 
ut  =  Table  [Evaluate  [Table  [M[i][t],  {i,  1,  md}]/.soln],  {t,  0,  tf}]; 
cc  =  Flatten[ut[[Dimensions[ut]  [[!]]]]]; 

Calculate  the  approximated  transformation. 

as[t_]  =  Sum[cc[[i]]p[t,f  —  1],  {«,  l,md}] 

-8.7142+1.09115t+0.146231(-2+2t2)-0.0114968(-6y2t+2y2t3)_Q_QQ^24268(12- 

24^2  +  At^) 

Calculate  the  actual  transformation  with  8192  samples. 

a  =  —6; 

6  =  6; 
n  =  13; 
nnd  =  2^3  ^  ij 


nd  =  Table[a  +  hk,  {k,  0,  nnd  —  1}]; 
odata  =  Table[0,  {«,  1,  nnd}]; 

odata  =  Flatten[Table[yl[a;]/.NDSolve[{yl'[t]  ==  y2[t], 

y2'|(]  ==  _(“±i^)((yl[t|)A3)  _ 

yl[0]  ==  10,y2[0]  ==  ■u},{yl,y2},{t,0,tf}]/.a;  ^tf,{i,l,nnd}]]; 
ListPlot[Transpose[{nd,  odata}],  PlotRange  ^  {—2,  —10}] 
Display ["nonlinsamptransactual.eps",  %,  "EPS"] 
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T  =  Transpose  [{nd,  odata}]; 

pp  =  Interpolation  [TjInterpolationOrder  1]; 

Plot  the  actual  and  approximate  transformations  and  calculate  the  weighted  error. 

Plot[{as[t], pp[t]},  {t,  —6, 6},  AxesLabel  — >  "g"}, 

PlotStyle  — >  {Dashing[{.03}],  GrayLevel[0]}] 

Eixport["gandghatiil4.eps",  %,  "EPS"] 

WeightedError  =  Sqrt[NIntegrate[(pp[t]  —  as[t])^/[t],  {t,  —6,6}]] 


0.0164006 

B.  2  Uniform  Sampling 

Dehne  the  parameters  and  conditions  of  the  nonlinear  ODE. 

<  <  Statistics‘  ContinuousDistributions‘ 

<  <  Graphics'  Graphics' 
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a  =  —6; 

6  =  6; 
kb  =  1.; 
kt  =  0.02; 
tf  =  10; 
m  =  10.; 
kl  =  10.; 
k3b  =  1; 
k3t  =  0.02; 

ti  =  0; 

npdf[t_]  =  PDF[NormalDistribution[0,  l],t]; 

/[t_]  =  npdf[t] 

e~^ 

Calculate  the  actual  transformation  with  8192  samples. 

n  =  13; 
nnd  =  2"  +  1; 
h  =  - 

nd  =  Tablefa  +  hk^  {fc,  0,  nnd  —  1}]; 
odata  =  Table[0,  {*,  1,  nnd}]; 

odata  =  Flatten[Table[yl[a;]/.NDSolve[{yl'[t]  ==  y2[t], 

y2'|«l  ==  -(“±i»“)((yl[t|)A3)  -  Byl[«|. 

yl[0]  ==  10,y2[0]  ==  u},{yl,y2},{t,0,tf}]/.a;  ^tf,{z,l,nnd}]] 

T  =  Transpose  [{nd,  odata}]; 

pp  =  Interpolation  [T,InterpolationOrder  1]; 

Plot[pp[t],  {t,  —6, 6},  AxesLabel  — >  {"x",  "g"}] 

Eixport[" nlactualtrans.eps",  %,  "EPS"] 
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g 


Sample  the  function  defining  the  transformation  from  ^  to  x(T,^)  .  Use  this  set  of 
samples  to  approximate  the  transformation. 

Using  8  Samples, 

n  =  3; 

mid  =  2”  +  1; 

nnd-l’ 

nd  =  Tablefa  +  hk,  {fc,  0,  nnd  —  1}]; 
odataS  =  Table[0,  {i,  1,  nnd}]; 

odataS  =  Flatten[Table[yl[x]/.NDSolve[{yl'[t]  ==  y2[t], 
y2'W  ==  -(“±if!D“)((yl[t]n3)  -  ttylW. 

yl[0]  ==  10,y2[0]  ==  ?;},{yl,y2},{t,0,tf}]/.a; ->tf,{z,l,nnd}]]; 

Interpolate  the  points  generated  from  the  uniform  sampling  method. 

T8  =  Transpose[{nd,  odataS}]; 

pp8  =  Interpolation[T8,InterpolationOrder  — »  1]; 

Plot[pp8[t],  {t,  -6,6}] 

Continue  above  steps  for  each  level  of  samples.  Plot  the  actual  and  approximate 
transformations  and  calculate  the  weighted  error  for  each  level  of  samples. 

Using  8  Samples, 
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Plot[{pp8[t],  pp[t]},  {t,  -6, 6}] 

WES  =  Sqrt[NIntegrate[(pp[t]  —  pp8[t])^/[t],  {t,  —6,6}]] 


0.127753 

Using  16  Samples, 

Plot[{ppl6[t],pp[t]},  {t,  — 6,6},  AxesLabel  ^  {"x",  "g”}, 
PlotStyle  — ^  {Dashing[{.03}],  GrayLevel[0]}] 

WE16  =  Sqrt[NIntegrate[(pp[t]  —  ppl6[t])^/[t],  {t,  —6,6}]] 


0.032066 

Repeat  this  process  for  the  rest  of  the  sample  levels.  Produce  the  error  plot  for  the 
nonlinear  ODE. 

LogLogListPlot[{{8,  WE8},  {16,  WE16},  {32,  WE32},  {64,  WE64},  {128,  WE128}, 

,  {256,  WE256},  {512,  WE512},  {1024,  WE1024},  {2048,  WE2048},  {4096,  WE4096}}, 
AxesLabel  ^  {"Number  of  Samples",  "Weighted  Error"},  Plot  Joined  — >  Ttue, 
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Ticks  — {{8, 32, 256, 1024, 4096},  Automatic}] 
Ebq)ort[" iilsamplingerror.eps",  %,  "EPS"] 


Weighted  Error 


B.  3  Nonuniform  Sampling 

Define  the  parameters  and  conditions  of  the  nonlinear  ODE. 

<  <  Statistics‘  ContinuousDistributions' 

<  <  Graphics'  Graphics' 

a  =  —6; 

6  =  6; 
kb  =  L; 
kt  =  0.02; 
tf  =  10; 
m  =  10.; 
kl  =  10.; 
k3b  =  1; 
k3t  =  0.02; 
u  =  0; 
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dist  =  NormalDistribution[0, 1]; 
pdf[x_]  =  PDF  [dist,  a;]; 
cdf[x_]  =  CDF[dist,x]; 

Calculate  the  actual  transformation  with  8192  samples. 

n  =  13; 
nnd  =  2”  +  1; 
h  =  - 

id  =  Tablefa  +  hk,  {fc,  0,  nnd  —  1}]; 
od  =  Table[0,  {i,  1,  nnd}]; 

od  =  Flatten[Table[yl[a;]/.NDSolve[{yl'[t]  ==  y2[f], 

ym  ==  -(!2t!aiilM‘)((yi[(|)A3)  _ 

yl[0]  ==  10,  y2[0]  ==  t;},  {yl, y2},  {t,  0,  tf}]/.x  ^  tf, 

{i,  l,nnd}]]; 

T  =  Transpose  [{id,  od}]; 

pp  =  Interpolation  [T,InterpolationOrder  1]; 


Calculate  the  inverse  cdf  of  the  normal  distribution. 

in  =  500; 

s  =  Table[iV[{cdf[-6  +  ^],  -6  +  ^}],  {i,  0,  in}]; 
ListPlot[s] 

cdfi  =  Interpolation[s,  InterpolationOrder  — >  1]; 
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Evaluate  the  inverse  cdf  to  get  a  nonuniform  sample  from  the  normal  distribution. 

Using  11  Samples, 

k  =  3; 
n  =  2^  —  1; 

XX  =  Tablef;^,  {i,  1,  n}]; 

2/  =  Table[0, 0,1,  n}]; 

y  =  cdfi[xx]; 

nd  =  Table [0,  {i,  l,n}]; 

nd  =  Flatten[{a,  MUM,  t/[[Length[t;])]+6^ 

{-6,  -3.57519,  -1.15037,  -0.674508, 

-  0.318658, 0.,  0.318658,  0.674508, 1.15037,  3.57519,  6} 

Sample  the  function  defining  the  transformation  from  ^  to  x(T,i^)  .  Use  this  set  of 
samples  to  approximate  the  transformation. 

nnd  =  Length  [nd] 

odata  =  Flatten[Table[yl[z]/.NDSolve[{yl'0]  ==  y20], 
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Til  =  Transpose[{nd,  odata}]; 

ppll  =  Interpolation[Tll,InterpolationOrder  — 1]; 

Plot[ppll[t],{t,  -6, 6}]; 

Continue  above  steps  for  each  level  of  samples.  Plot  the  actual  and  approximate 
transformations  and  calculate  the  weighted  error  for  each  level  of  samples. 

Using  11  Samples, 

Plot[{ppll[t],pp[t]},  {t,  -6,6}, 

PlotStyle  — >  {Dashing[{.03}],  GrayLevel[0]}] 

WEll  =  Sqrt[NIntegrate[(pp[t]  —  ppll[t])^pdf[t],  {t,  —6,6}]] 


0.130374 

Using  19  Samples, 

Plot[{ppl9[t], pp[t]},  {t,  -6, 6}, 

PlotStyle  {Dashing[{.03}],  GrayLevel[0]}] 

WE19  =  Sqrt[NIntegrate[(pp[t]  -  ppl9[t])2pdf[t],  {t,  -6,6}]] 
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0.0722684 


Repeat  this  process  for  the  rest  of  the  sample  levels.  Produce  the  error  plot  for  the 
nonlinear  ODE. 

LogLogListPlot[{{ll,  WEll},  {19,  WE19},  {35,  WE35},  {67,  WE67}, 

{131,  WE131},  {259,  WE259},  {515,  WE515},  {1027,  WE1027}, 

{2051,  WE2051},  {4099,  WE4099}}, 

AxesLabel  ^  {"Number  of  Samples",  "Weighted  Error"},  Plot  Joined  — True, 

Ticks  {{11, 35, 259, 1027, 4099},  Automatic}] 

Eixport[" iilnonuiiiformsamplingerror.eps",  %,  "EPS"] 


Weighted  Error 
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Appendix  C.  Mathematica  Code  for  the  CDF 
Transformation  Method 


Define  the  parameters  and  conditions  of  the  linear  ODE. 

<  <  Statistics‘  ContinuousDistributioiis‘ 

<  <  Graphics'  Graphics' 

npdf[t_]  =  PDF[NormalDistribution[0, 

/[t_]  =  npdf[t]; 
a  =  —6; 

6  =  6; 
pd  =  4; 
tf  =  10; 
m  =  10; 
kb  =  10; 
kt  =  2; 
c  =  4; 

V  =  25; 

Create  A  and  B  matrices  for  S. 


ipa  =  IdentityMatrix[pd  +  1]; 

MatrixForm[ipa] ; 

ipb  =  Table[0,  {i,  0,  pd},  {j,  0,  pd}]; 

For[z  =  0,  *  <  pd,  i++, 

For[7  =  i,j  <  pd,  j++, 

ipb[[«  +  l,i  +  1]]  =  Chop[NIntegrate[tp[t,*]p[t,  j]/[t],  {t,  -Infinity,  Infinity}]] 
vph[\j  +  l,i  +  1]]  =  ipb[[i  +  l,i  +  1]]; 

] 

] 
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MatrixForm  [ipb] ; 

md  =  Dimensions[ipa][[l]]; 

aall  =  Table[0,  {i,  1,  md},  {j,  1,  md}]; 
aal2  =  IdentityMatrfac[md]; 

aa21  =  — ((^)IdentityMatrix[md]  +  (^)Inverse[ipa].ipb); 
aa22  =  —  ^IdentityMatrix[md]; 

aa  =  Chop[Join[Transpose[Join[aall,  aal2]],Transpose[Join[aa21,aa22]]]] 
MatrbcForm[aa] ; 

Solve  the  system  of  linear  ODEs. 

ua  =  Table  [M[i][f],  {i,  1, 2md}]; 
upa  =  Table[M[i]'[t],  {i,  l,2md}]; 

ics  =  Flatten  [{Table  [w[i][0]  ==  0,  {i,  1,  md}],  wfmd  +  1][0]  ==  v, 
Table[w[i][0]  ==  0,  {i,  md  +  2, 2md}]}]; 
eqns  =  Flatten[Table[upa[[z]]  ==  (aa.ua)[[i]],  {*,  1, 2md}]]; 
eqns  =  Flatten[{eqns,ics}]; 

MatrbcForm[eqns] ; 

soln  =  NDSolve[eqns,  Tablefttfi],  {i,  1, 2md}],  {t,  0,  tf}]; 
ut  =  Table  [Evaluate  [Table  [M[i][t],  {i,  1,  md}]/.soln],  {t,  0,  tf}]; 

Table  [Evaluate  [Table  [M[i][t],  {f,md  +  1, 2md}]/.soln],  {t,  0,tf}]; 
cc  =  Flatten[ut[[Dimensions[ut]  [[!]]]]] 

(-0.537323,  -2.02877,  0.514555,  0.941688,  -0.584251} 

Calculate  the  approximated  transformation. 

as[t_]  =  Sum[cc[[i]]p[t,i  —  1],  {i,  l,md}]; 

Plot[as[t],{t,  — 6,6},AxesLabel  — >  {"x",  "ghat(x)"}] 
nd  =  2^2  +  1; 

id  =  Table[a  +  (i  -  1)^,  {*,  1,  nd}]; 
od  =  Table[0,  {i,  1,  nd}]; 

For[i  =  1,  i  <  nd,  Z++, 
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od[[z]]  =  as[id[[*]]] 

] 


Using  approximate  transformation  just  found,  generate  the  pdf  using  the  transforma¬ 
tion  method. 

exlim  =  Table  [{0, 0, 0, 0, 0},  {i,  1,  nd}]; 
cnt  =  1; 

exlim[[cnt]]  =  {1,  id[[cnt]],od[[cnt]],  1,0}; 
index  =  2; 

If[od[[index]]  —  od[[index  —  1]]  <  0,  dir  =  —1,  dir  =  1]; 

While  [index  <  nd, 
index-h-l-; 

S  =  Sign[od[[index]]  —  od[[index  —  1]]]; 

If[(5  ==  dir, , 
cnt-l— 1-; 

exlim[[cnt]]  =  {index  —  1,  id[[index  —  1]],  od[[index  —  1]],  1, 0}; 
dir  =  — Sign[dir] 

]; 

] 

cnt-|— 1-; 

exlim[[cnt]]  =  {index,  id[[index]],od[[mdex]],  1,0}; 
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exlim  =  Take[exlim,  cnt]; 

If[exlim[[l,  3]]  >  exlim[[2,3]], 
exlim[[l,  5]]  =  1;  exlim[[2, 5]]  =  —1, 
exlim[[2, 5]]  =  1;  exlim[[l,  5]]  =  —1 

]; 

For[i  =  3,  *  <  cnt,  i++, 

exlim[[i,  5]]  =  If[exlim[[i  —  1, 5]]  ==  1,  —1, 1] 

]; 

nex  =  Dimensions[exlim]  [[!]]; 

If[od[[index]]  —  od[[index  —  1]]  <  0, 
max  =  Take[exlim,  {1,  cnt,  2}]; 
min  =  Take[exlim,  {2,  cnt  —  1, 2}], 
min  =  Take[exlim,  {1,  cnt,  2}]; 
max  =  Take[exlim,  {2,  cnt  —  1, 2}] 

]; 

max  =  Sort  [max,  #1  [[2]]  <  #2[[2]]&:]; 
min  =  Sort[min,#l[[2]]  <  #2[[2]]&]; 
nmax  =  Dimensions[max]  [[!]]; 
nmin  =  Dimensionsfmin]  [[!]]; 

node  =  iV[Transpose[{Table[i,  {«,  1,  nd}],  id,  od,  Table[0,  {*,  1,  nd}]}]] 
ysnode  =  Sort[node,  #1[[3]]  <  #2[[3]]&:]; 
exlim  =  Chop  [exlim] 

(^remove  repeated  yvalnes  from  ysnode*) 
new  =  {ysnode[[l]]}; 

For[i  =  l,i  <  Dimensions[ysnode] [[!]], i++, 

If[Chop[ysnode[[i,  3]]==ysnode[[i  +  1, 3]]], 

(*do  nothing*), 

(*add  to  new*)new  =  Join[new,  {ysnode[[i  +  1]]}] 
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] 

] 

{{1,  -6,  -180.908, 1,  -1},  {1408,  4.49415, 1, 1}, 

(2430,  -3.11646, 1,  -1},  (3135,  -0.295545, 1, 1},  {4097,  6,  -53.0141, 1,  -1}} 

For  [(*  start  at  i  =  2  because  i  =  1  is  the  miniTniim  in  new  so  there  is  no  contribution  , 
then  process  each  element  of  new  *) 
i  =  2,i  <  Dimensions[new][[l]],i++, 

j  =  1; 

(*nex  is  the  number  of  entries  in  exlim  ,  which  contains  the  local  extrema 
and/or  the  limits  of  integration  *) 

While  \j  <  nex, 

(*  Print  ["j  =  ",j,  "  nex  =  ",nex];*) 

If[(*  an  extreme  *)exlim[[j,  4]]  ==  1, 

If[(*a  minimum*)exlim[[j,  5]]  ==  —1, 

(*Print["i  =  ",z,  "  j  =  ",  j,  "  min"];*) 

If[new[[*,  3]]  >  exhm[[7,3]], 

(*add  limits  and  remove  min  *) 
tm  =  exlim[[)]]; 

If[tm[[l]]  i-  1, 

xtmp  =  (node[[tm[[l]]  —  1,2]]  —  node[[tm[[l]],2]])/(node[[tm[[l]]  —  1,3]]  —  node[[tm[[l]],3]]) 
(new[[i,  3]]  -  node[[tm[[l]],  3]])  +  node[[tm[[l]],  2]]; 
tmp  =  {{tm[[l]],xtmp,  new[[z,3]],0,  — 1}}; 

exlim  =  Flatten[{Take[exhm,  J  —  l],tmp,  Take[exlim,  j  —  nex]},  1], 

tmp  =  {{tm[[l]],node[[tm[[l]],2]],new[[i,3]],0,-l}}; 

exhm  =  Flatten[{Take[exhm,j  —  l],tmp,Take[exlim,  j  —  nex]},  1] 

]; 

(*Print  [exlim];  *) 

If[tm[[l]]  ^  nd, 
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xtmp  =  (node[[tm[[l]]  +  1,2]]  —  node[[tm[[l]],  2]])/(node[[tm[[l]]  +  1,3]]  —  node[[tm[[l]],  3]]) 
(new[[i,  3]]  -  node[[tm[[l]],  3]])  +  node[[tm[[l]],  2]]; 
tmp  =  {{tm[[l]],xtmp,new[[«,3]],0, 1}}; 

exlim  =  Flatteii[{Take[exlim,j],  tmp.  Take  [exlim,j  —  nex]},  1], 

tmp  =  {{tm[[l]],node[[tm[[l]],2]],new[[i,3]],0, 1}}; 

exlim  =  Flatten[{Take[exlim,  j],  tmp.  Take  [exlim,  j  —  nex]},  1] 

]; 

(*Print  [exlim];  *) 
nex++ 

* 

(*do  nothing*) 

] 

,  (*else  a  maximum*) 

(*If[(i  >  192)&:&:(i  <  195),Print["i  =  ''  j  =  '',j,  "  max",  "  nex  =  ",nex]; 

Print  [exlim]; 

Print[new[[i,  3]],  "  ", exlim[[j,  3]]], 

];*) 

If[new[[i,  3]]  >  exlim[[7,3]], 

(*removelimitsand  max  *) 

If[(nex  >  3)&:&:(j  l)&:&(j  ^  nex), 

exlim  =  Flatten[{Take  [exlim,  j  —  2],  Take[exlim,  j  +  1  —  nex]},  1]; 
nex  =  nex  —  3; 
j=j-  3, 

If[(i  ^  l)&&(i  ^  nex), 
exlim  =  Take  [exlim,  —2]; 
nex  =  2, 

ifb  ==  1, 

exlim  =  Flatten[{Take  [exlim,  1],  Take[exlim,  2  —  nex]},  1]; 
nex  =  nex  —  1; 
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exlim[[l,4]]  =  0; 
exlim[[l,  5]]  =  —  1, 

If|j  ==  nex, 

exlim  =  Flatten[{Take[exlim,  nex  —  2],  Take[exlim,  —1]},  1]; 
nex  =  nex  —  1; 
exlim[[nex,  4]]  =  0; 
exlim[[nex,  5]]  =  1, 

](*endlf[7  ==  nex...*) 

](*endlf[7  ==  1...*) 

](*Ifb-  ^  1...*) 

](*endlf[nex  >  3...*) 

(*else  do  nothing*) 

](*end  If[new[[i,  3]]...*) 

](*end  If  (a  minimum  )*) 

];  (*end  If  (an  extreme  )*) 

j++ 

];  (*end  Whileb  <  nex*) 

(*Interpolate  to  define  limits  of  integration*) 

Forb  =  l,i  <  nex,  j++, 

If[(*a  limit*)exlim[b,  4]]  ==  0, 

If[(*left  limit*)exlim[b,  5]]  ==  —1, 

While[((exlim[b,  1]]  >  l)&:&:(new[[i,  3]]  >  node[[exlim[b,  1]]  —  1,3]])), 
exlim[b,  1]]  =  exlim[b,  1]]  -  1 
](*end  While[((exlim...*); 
in  =  exlim[b,  1]]; 

If[in  ^  1, 

exlimfb,  2]]  =  (node[[in  —  1, 2]]  —  node[[in,  2]])/(node[[in  —  1, 3]]  —  node[[in,  3]]) 
(new[[i,  3]]  —  node[[in,  3]])  +  node[[in,  2]]; 
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exlim[[j,  3]]  =  new[[i,3]], 
exlim[[j,  2]]  =  node[[in,  2]]; 
exlim[[j,  3]]  =  new[[*,  3]] 

](*end  If  [in  ^  1*) 

,  (*else  a  right  limit*) 

While[((exlim[[7, 1]]  <  nd)&&(new[[i,  3]]  >  node[[exlim[[j,  1]]  +  1, 3]])), 
exlim[[j,  1]]  =  exlim[[7, 1]]  +  1 
]; 

in  =  exlim[[j,  1]]; 

If[in  ^  nd, 

exlim[[j,  2]]  =  (node[[in  +  1, 2]]  —  node[[in,  2]])/(node[[in  +  1, 3]]  —  node[[in,  3]]) 

(new[[i,  3]]  —  node[[in,  3]])  +  node[[in,  2]]; 

exlim[[j,3]]  =  new[[i,3]], 

exlim[[j,2]]  =  node[[in,2]]; 

exlim[[j,  3]]  =  new[[i,  3]] 

] 

](*end  If[(*left  limit*)...*) 

(*do  nothing*) 

](*end  If[(*a  Umit*)...*) 

];  (*end  For\j  =  1...*) 

(*integrate  and  update  cdf*) 

j  =  1; 

While  [)  <  nex, 

While[((exlim[[7,4]]  ^  0)&:&(j  <  nex)), 

j++ 

]; 

(*If[i  >  220,  Print  [''integrate  ",  "i  =  ",i,  "  j  =  ",i,  "  ",new[[*,4]],  "  ",exlim],];  *) 
If[exlim[[j,4]]  ==  0, 
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new[[i,  4]]  =  new[[z,  4]]  +  NIntegrate[npdf[a;],  {x, exlim[[_7, 2]],  exlim[|j  +  1, 2]]}], 

]; 

i  =  i  +  2 

] 

]; 

(*end  For  i*) 

cdf  =  Transpose  [Take  [Transpose  [new],  —2]]; 

ListPlot[cdf] 


Differentiate  the  cdf  to  get  the  pdf. 

pdf  =  Table[{cdf[[*,  1]],  (cdf[[z  +  1, 2]]  -  cdf[[z  -  1, 2]])/(cdf[[*  +  1, 1]]  -  cdf[[z  -  1, 1]])  }, 
{*,  2,  Dimensions[cdf][[l]]  —  1}]; 

pp  =  Interpolation  [pdf,  InterpolationOrder  — >  1]; 

Plot[{pp[t]},{t,  -6, 6}] 
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This  is  the  actual  transformation  for  the  model  problem. 

qqls-t-]  =  -  4(^)(kb + skt)]]; 

=  sqrt[(^)^-44)(kb+.kt)]^^P[-^]Skih[fSqrt[(^)"  -  4(^)(kb  +  skt)]]; 

Plot[{g[t],  as[t]},  {t,  —6, 6},  AxesLabel  {"x",  "g"},  PlotStyle  —*■  { Gr  ay  Level  [0],  Dashing[{.03}]}] 

id  =  Table[a  +  (i  -  1)^,  {i,  1,  nd}]; 
od  =  Table[0,  {i,  1,  nd}]; 

For[i  =  1,*  <  nd,  Z++, 
od[W]  =  g[id[[*]]] 

] 


g 


Find  the  cdf  of  the  actual  transformation  with  the  transformation  method  using  the 
same  code  used  to  approximate  the  transformation. 
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Plot  the  actual  and  approximate  pdfs  and  calculate  the  error. 

Plot[{pp[t],ppl[t]},  {t,  -6,6},  AxesLabel  {"y",  "/y"}, 
PlotStyle  — >  {GrayLevel[0],Dashing[{.03}]}] 

Display ["lingaiidghat4pdf.eps",  %,  "EPS"] 

WeightedError  =  Sqrt[NIntegrate[(pp[t]  —  ppl[t])^,  {t,  —3,4.1}]] 

fy 


0.0228001 
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Appendix  D.  Mathematica  Code  for  Monte  Carlo  Simulation 


Define  the  parameters  and  conditions  of  the  linear  ODE. 

<  <  Statistics‘  ContinuousDistributioiis‘ 

<  <  Graphics'  Graphics' 

npdf[t_]  =  PDF[NormalDistribution[0, 

/[t_]  =  npdf[t]; 
a  =  —6; 

6  =  6; 
pd  =  4; 
tf  =  10; 
m  =  10; 
kb  =  10; 
kt  =  2; 
c  =  4; 

V  =  25; 

Create  A  and  B  matrices  for  S. 


ipa  =  IdentityMatrix[pd  +  1]; 

MatrixForm[ipa] ; 

ipb  =  Table[0,  {i,  0,  pd},  {j,  0,  pd}]; 

For[z  =  0,  *  <  pd,  i++, 

For[7  =  i,j  <  pd,  j++, 

ipb[[«  +  l,i  +  1]]  =  Chop[NIntegrate[tp[t,«]p[t,  j]/[t],  {t,  -Infinity,  Infinity}]]; 
vph[\j  +  l,i  +  1]]  =  ipb[[«  +  l,i  +  1]]; 

] 

] 

MatrixFor  m  [ipb] ; 

md  =  Dimensions[ipa][[l]]; 
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aall  =  Table[0,  {z,  1,  md},  {j,  1,  md}]; 
aal2  =  IdentityMatrix[md]; 

aa21  =  — ((^)IdentityMatrix[md]  +  (^)Inverse[ipa].ipb); 
aa22  =  —  ^IdentityMatrix[md]; 

aa  =  Chop[Join[Traiispose[Join[aall,  aal2]],Traiispose[Join[aa21,aa22]]]]; 
Mat  rixForm  [aa] ; 

Solve  the  system  of  linear  ODEs. 

ua  =  Table  [M[i][t],  {i,  1, 2md}]; 
upa  =  Table[M[z]'[f],  {i,  l,2md}]; 

ics  =  Flatten[{Table[M[i][0]  ==  0,  {i,  1,  md}],  wfmd  +  1][0]  ==  v, 
Table[M[z][0]  ==  0,  {i,  md  +  2, 2md}]}]; 
eqns  =  Flatten[Table[upa[[«]]  ==  (aa.ua)  [[«]],  {«,  1, 2md}]]; 
eqns  =  Flatten[{eqns,ics}]; 

MatrixForm[eqiis] ; 

soln  =  NDSolvefeqns,  Table[u[i],  {i,  1, 2md}],  {t,  0,  tf}]; 
ut  =  Table  [Evaluate  [Table  [M[i][t],  {z,  1,  md}]/.soln],  {t,  0,  tf}]; 

Table  [Evaluate  [Table  [M[i][t],  {z,md  +  1, 2md}]/.soln],  {t,  0,tf}]; 
cc  =  Flatten[ut[[Dimensions[ut]  [[!]]]]] 

{-0.537323,  -2.02877,  0.514555,  0.941688,  -0.584251} 

Calculate  the  approximated  transformation. 

as[t_]  =  Sum[cc[[i]]j5[t,f  —  1],  {*,  l,md}]; 

Plot[as[t],{t,  — 6,6},AxesLabel  {"x",  "ghat(x)"}] 
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ghat (x) 


Using  approximate  transformation  just  found,  generate  the  pdf. 

ss  =  30000; 

ri  =  Sort[Table[Random[NormalDistribution[0, 1]], 
data  =  Table  [0,  {i,  l,ss}]; 

For[i  =  l,i  <  ss,i++, 
data[[i]]  =  as[ri[[i]]] 

] 

aa  =  Histogram  [data,  HistogramScale  — » 1,  AxesLabel  — >  {"y",  "fy"}] 


fy 


Create  a  pdf  function  by  interpolating  the  midpoints  of  the  top  of  each  rectangle  of 
the  histogram. 

aaa  =  Transpose[Transpose[aa[[l]]]  [[!]]]  [[2]]; 
datal  =  Table[{0, 0},  {i,  1,  Length[aaa]  +  2}]; 
datal[[l]]  =  {aaa[[l,  1,  l]],aaa[[l,2,2]]/2}; 

For|j  =  l,i  <  Lengthfaaa],  j++, 
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U  =  aaa[[7, 1]]; 
ur  =  aaa[[7,2]]; 

datal|b  +  lll  =  {™l¥®ll.ur[[2]]} 

1 

datal[[Length[datal]]]  =  {aaa[[Length[aaa],2]][[l]],0}; 
ah  =  Interpolation[datal,  InterpolationOrder  1]; 
Plot[{ah[a;]},  {;r,  —5, 5}] 


This  is  the  actual  transformation  for  the  model  problem. 

qq[s-.t-]  =  -  4(^)(kb  +  skt)]]; 

9[s-]  =  sqrt[(^)2-44)(kb+akt)]^^P[-^]^^^^[fSqrt[(^)"  -  4(i)(kb  +  skt)]]; 
nd  =  2®  +  1; 

Plot[{q'[f],  as[f]},  {t,  —6, 6},  AxesLabel  — >  {"x",  "g"},  PlotStyle  — »  {GrayLevel[0],  Dashing[{.03}]}] 
id  =  Table[a  +  (i  -  1)^,  {i,  1,  nd}]; 
od  =  Table[0,  {i,  1,  nd}]; 

For[«  =  1,«  <  nd,«++, 
od[[i]]  =  g[id[[i]]] 

] 
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g 


Find  the  pdf  of  the  actual  transformation  using  the  transformation  method  as  de¬ 
scribed  in  the  previous  appendix.  Plot  the  actual  and  approximate  pdfs  and  calculate 
the  error. 

Plot[{ah[t],ppl[t]},  {t,  — 6,6},  AxesLabel  — {"y",  "fy"}, 

PlotStyle  ^  {GrayLevel[0],Dashing[{.03}]}] 

Display["MChermite4pdf.eps'',  %,  "EPS"] 

WeightedError  =  Sqrt[NIntegrate[(ah[t]  —  ppl[f])^,  {t,  —3,4.1}]] 

fy 


0.0323146 
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